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Abstract 

We  study  a  finite-difference  discretization  of  an  ill-posed  nonlinear  parabolic  par¬ 
tial  differential  equation.  The  PDE  is  the  one-dimensional  version  of  a  simplified 
two-dimensional  model  for  the  formation  of  shear  bands  via  anti-plane  shear  of  a 
granular  medium.  For  the  discretized  initial  value  problem,  we  derive  analytically, 
and  observed  numerically,  a  two-stage  evolution  leading  to  a  steady-state:  (i)  an 
initial  growth  of  grid-scale  instabilities,  and  (ii)  coarsening  dynamics.  Elaborating 
the  second  phase,  at  any  fixed  time  the  solution  has  a  piecewise  linear  profile  with 
a  finite  number  of  shear  bands.  In  this  coarsening  phase,  one  shear  band  after  an¬ 
other  collapses  until  a  steady-state  with  just  one  jump  discontinuity  is  achieved. 
The  amplitude  of  this  steady-state  shear  band  is  derived  analytically,  but  due  to 
the  ill-posedness  of  the  underlying  problem,  its  position  exhibits  sensitive  depen¬ 
dence.  Analyzing  data  from  the  simulations,  we  observe  that  the  number  of  shear 
bands  at  time  t  decays  like  t  From  this  scaling  law  we  show  that  the  time-scale 
of  the  coarsening  phase  in  the  evolution  of  this  model  for  granular  media  critically 
depends  on  the  discreteness  of  the  model. 

Our  analysis  also  has  implications  to  related  ill-posed  nonlinear  PDEs  for  the  one¬ 
dimensional  Perona-Malik  equation  in  image  processing  and  to  models  for  clustering 
instabilities  in  granular  materials. 
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1  Introduction 


The  degenerate  parabolic  PDE, 


(1.1) 


in  which  R  is  a  rotation  matrix,  arises  from  a  simplified  model  of  the  velocity  field  of  a 
sheared  granular  material  [33].  This  equation  is  ill-posed,  a  property  typical  of  continuum 
models  of  granular  media  [25,31,32,34],  To  study  the  dynamics  of  this  model,  in  this  paper 
we  analyze  finite  difference  approximations  of  the  PDE  (1.1).  More  precisely,  we  discretize 
in  space,  and  study  the  resulting  system  of  ODEs,  which  we  refer  to  as  the  discrete  model. 
This  model  is  a  form  of  regularization  of  (1.1);  the  discrete  model  is  well-posed,  mollifying 
instabilities  at  the  highest  frequencies. 


The  simulations  show  that  amplification  of  perturbations  at  short  wavelengths  lead  quickly 
to  small-amplitude  discontinuities,  see  Figure  lb.  This  short  time  behavior  is  followed  by 
a  gradual  coarsening  and  long  time  evolution  to  a  steady  state.  The  steady  state  has  a 
single  discontinuity  that  we  refer  to  as  a  shear  band  [9],  see  Figure  lc.  The  time-scale  for 
the  coarsening  dynamics  phase  of  the  evolution  diverges  as  the  mesh  scale  Ax  tends  to 
zero.  Thus,  strictly  speaking,  in  the  continuum  limit,  no  evolution  occurs!  In  other  words, 
the  discreteness  of  the  finite  difference  model  is  essential  to  the  dynamics  of  this  system. 
Implications  of  this  result  for  granular  flows  will  be  discussed  further  in  a  subsequent  paper. 


The  behavior  of  this  nonlinear  ill-posed  equation  should  be  contrasted  with  the  behavior  of 
ill-posed  linear  equations  where  short  wavelength  disturbances  are  simply  amplified  catas¬ 
trophically.  Correspondingly,  solutions  of  difference  approximations  of  linear  ill-posed  equa¬ 
tions  diverge  in  every  norm,  as  the  mesh  spacing  approaches  zero.  For  our  discrete  model 
of  (1.1),  as  Ax  — »  0,  the  sequence  of  solutions  of  the  difference  equations  converges  in  L'2, 
but  because  of  ill-posedness  it  blows-up  in  H1.  The  divergence  in  H 1  is  a  consequence  of  the 
very  rapid  development  of  “infinitesimal  discontinuities”  in  the  solution.  These  observations 
on  the  dynamics  in  the  system  are  fully  developed  in  Section  6. 


The  linear  ill-posedness  of  the  time-dependent  equations  is  related  to  the  property  that 
the  steady-state  equations  are  hyperbolic.  The  steady-state  equations  are  well-known  to 
support  shock  waves  [22],  In  nonlinear  continuum  descriptions  of  granular  materials  linear 
ill-posedness  means  that  initial  value  problems  are  sensitive  to  small  perturbations,  but  it 
also  provides  a  mechanism  for  the  formation  of  fine-scale  localized  structures  such  as  shear 
bands  and  shocks. 


PDE  models  for  granular  flow  are  the  result  of  treating  the  material  as  a  continuum,  allowing 
dynamics  at  all  wavelengths,  whereas  in  fact  wavelengths  much  shorter  than  the  grain  size  of 
particles  in  the  granular  media  have  no  physical  significance.  Regarding  the  finite  difference 
mesh  parameter  as  being  on  the  order  of  the  grain  size,  we  effectively  replace  the  continuum 
description  by  a  discrete  model  that  may  more  faithfully  represent  the  range  of  wavelengths 
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Fig.  1.  One-dimensional  evolution  in  (1.1);  (a)  initial  data,  (b)  short  time  behavior,  (c)  long-time 
evolution  to  the  steady-state  shear  band. 

relevant  for  granular  flow.  The  discrete  model  is  analogous  to  discrete  mechanical  models 
whose  continuum  limit  yields  a  continuum  description,  as  has  been  used  in  constructing 
continuous  models  of  high  density  discrete  mechanical  systems  [26-28]. 

Ill-posed  nonlinear  parabolic  equations  also  arise  in  population  dynamics  in  mathematical 
biology  [21,23],  edge  enhancement  in  image  processing  [2,5,8,20],  and  other  problems  in 
granular  media  flows  [10,19,35].  Since  the  early  work  of  Perona  and  Malik  [24],  ill-posed 
nonlinear  diffusion  equations  have  been  used  to  produce  enhancement  of  edges  in  digitized 
images  by  selective  amplification  of  intensity  gradients  via  backward  diffusion.  These  models 
have  attracted  a  lot  of  attention  in  the  engineering  and  mathematical  literature  [20],  and 
there  are  still  many  mathematical  open  questions.  While  (1.1)  has  a  very  different  motivation 
than  the  Perona-Malik  models,  it  shares  many  structural  similarities,  and  some  of  our  results 
may  have  further  implications  for  image  processing.  Recent  studies  [10,35]  of  clustering 
instabilities  (also  called  inelastic  collapse)  in  dilute  granular  gases  [19]  have  yielded  models 
with  non-monotone  flux  functions  and  discrete  diffusive  coupling  in  space.  The  discrete  model 
considered  here  has  some  qualitative  features  similar  to  the  clustering  dynamics  observed  in 
those  papers. 


1.1  Outline  of  paper 


The  remainder  of  this  paper  is  organized  as  follows.  In  Section  2  we  formulate  the  PDE 
model  in  one  and  two  space  dimensions,  and  give  preliminary  information  concerning  the 
discrete  model.  Specifically,  in  Section  2.1,  we  discuss  the  PDE  in  the  context  of  anti-plane 
shear,  and  specify  the  nature  of  the  linear  ill-posedness.  In  Section  2.2,  we  consider  plane 
wave  perturbations  of  uniform  shear.  We  show  that  the  resulting  PDE  in  one  space  variable 
and  time  has  the  character  of  a  backward-forward  heat  equation.  In  Section  2.3,  we  write 
down  the  discrete  model,  and  show  that  solutions  remain  bounded  globally  in  time,  and 
uniformly  in  mesh  spacing  Ax. 

In  Section  3,  we  describe  all  equilibrium  solutions  of  the  discrete  model.  Apart  from  the 
trivial  solution,  corresponding  to  uniform  shearing,  equilibrium  solutions  have  one  or  more 
finite  jumps  (that  persist  under  mesh  refinement),  which  we  call  shocks  or  shear  bands.  In 
Section  4,  we  show  that  the  trivial  solution  and  single  shear  band  solutions  are  the  only 
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possible  stable  equilibria.  We  characterize  precise  conditions  under  which  these  equilibria 
are  in  fact  stable.  Multiple  shock  solutions  are  shown  to  be  unstable,  a  property  that  helps 
explain  the  coarsening  exhibited  in  Figure  1.  The  proof  is  based  on  identifying  a  Liapunov 
function  for  the  discrete  model,  and  exploring  the  nature  of  it’s  critical  points. 

In  Section  5,  we  describe  detailed  dynamic  simulations  like  those  of  Figure  1.  In  Section  6, 
we  discuss  the  continuum  limit  Ax  — >-  0,  demonstrating  the  existence,  via  numerical  experi¬ 
ments,  of  a  scaling  law  for  the  evolution  of  shear  bands  in  the  coarsening  process.  Finally,  in 
Section  7,  we  analyze  the  coarsening  that  occurs  in  the  intermediate  dynamics.  Specifically, 
in  Section  7.1,  we  analyze  the  coarsening  from  a  solution  with  K  shocks  to  a  solution  with 
K  —  1  shocks,  when  K  is  large.  In  Section  7.2,  we  formulate  a  reduced  discrete  model  that 
describes  the  evolution  of  the  shocks  in  isolation  from  the  smooth  part  of  the  solution.  The 
comparison  of  predictions  from  the  reduced  model  and  the  full  model  helps  to  justify  our 
explanation  of  how  the  coarsening  process  takes  place. 


2  Formulation  of  the  problem 


2.1  The  continuous  two-dimensional  model 


The  dependent  variable  v  =  v(x,y,t )  in  equation  (1.1)  may  be  thought  of  as  the  velocity 
in  the  ^-direction  in  a  block  of  material  undergoing  anti-plane  shearing.  The  equation  for 
conservation  of  momentum  equates  acceleration,  dtv,  at  constant  density  (normalized  to 
unity)  with  the  divergence  of  stress,  V  •  r.  Modeling  stress  by  the  vector  t  =  RVv/\ V«j,  in 
which  R  is  the  matrix  representing  a  rotation  counterclockwise  through  a  constant  angle  a, 


R  = 


cos  a  -sma 
sin  a  cos  a 


(2.1) 


with  0  <  a  <  7r/2,  we  arrive  at  equation  (1.1).  This  model  and  related  equations  are  studied 
in  a  series  of  papers  [3,4,6,13,12,25,31,32,34],  Further  details  relating  this  constitutive  relation 
to  the  full  stress  tensor  and  to  properties  of  a  perfectly  plastic  material  are  given  in  [33]. 
Equation  (1.1)  is  a  special  case  of  the  model  in  [33]  subject  to  constant  yield  strength, 
specifically  |r|  =  1. 


The  following  proposition,  adapted  from  [33],  identifies  the  ill-posedness  exhibited  by  (1.1). 

Consider  a  specific  solution  v® (x,y,t)  of  (1.1).  If  we  linearize  the  equation  about  vq  and 

freeze  coefficients  of  the  linearized  equation  at  a  point  (xo,yo,to),  the  resulting  equation 

— # 

admits  solutions  in  the  form  exp{At  +  i(£\x  +  £2 y)}-  The  equation  relating  A  to  £  =  (£1,^2) 

is  the  dispersion  relation  for  the  linear  equation.  For  this  equation,  A  is  real,  and  we  say  the 

— > 

equation  is  ill-posed  if  A  is  not  uniformly  bounded  above  in  £.  The  proposition  describes  the 

— # 

nature  of  this  ill-posedness,  identifying  a  wedge  of  directions  £  within  which  A  is  positive 
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— # 

and  unbounded,  approaching  infinity  quadratically  in  £. 

Proposition  1  A  solution  v0(x,  y,  t )  of  (1.1)  is  linearly  ill-posed  in  a  neighborhood  of(x 0,  yo,  t0) 

— #  — > 

with  respect  to  the  exponential  perturbation  exp^^x  +  £2 y)}  if  £  (or  —£)  lies  in  the  sector 

arg(Vno)  <  arg(|)  <  arg(Vu0)  +  a  (2.2) 


where 


arg(Vn0)  =  arctan 


(  dyVo(x0,y0,to) 

\dxv0(x0,y0,t0) 


Note  that  the  upper  bound  in  (2.2)  is  the  argument  of  the  rotated  gradient  vector,  arg(i?Vn0)  = 
arg(Vuo)  +  a. 


Proof:  Using  the  summation  convention,  we  write  out  (1.1)  as 


dv 

dt 


1 

Vn 


djdkv  - 


dhff 

Vu 


(2.3) 


Substitute 

v  =  v0(x,  y,t)  +  e  et^lX+^eM 

into  (2.3),  equate  terms  of  order  e,  and  freeze  coefficients.  Then,  writing 

MO  =  \(0  +  o(\£\) 

to  extract  the  principal  part  in  the  growth  rate 1 ,  we  calculate  that 

Ap(f )  =  ||f  cos  a  -  ( Rd ,  0(«,  £ )  =  (a  x  £ )  •  (Rd  x  £ )  (2.4) 

where  (•,  •)  denotes  the  Euclidean  inner  product  and  a  =  Vn0/|Vn0|.  The  second  equality  in 

(2.4)  is  a  consequence  of  Lagrange’s  vector  identity  and  yields  the  growth  rate  in  terms  of 

— > 

vector  dot  and  cross  products.  Observe  that  Xp  vanishes  if  £  is  parallel  to  a  or  Ra.  A  nonzero, 

homogeneous  quadratic  form  such  as  (2.4)  can  vanish  along  only  two  directions.  Therefore 

A p(£)  cannot  change  sign  within  the  wedge  (2.2).  We  may  see  that  Ap(£)  is  positive  in  this 

— > 

wedge,  and  negative  outside  it,  by  observing  that  on  the  circle  {|£|  =  1},  the  first  term  in 

(2.4)  is  constant  while  the  second  term  assumes  its  maximum  where  the  circle  intersects  the 

— # 

line  that  bisects  the  wedge,  {arg£  =  arga  +  ct/2}.  ■ 

As  may  be  seen  from  the  proof,  the  two  angles  that  bound  the  wedge  (2.2)  represent  char¬ 
acteristic  directions  in  the  (x,  y)-plane  of  the  steady-state  equation 

div  (fiiw) =  0;  <2'5) 

1  Note  that  if  vq (x,y)  is  an  equilibrium  solution  of  (1.1),  then  A  =  \p.  i.e.  the  exponential  growth 
rate  equals  the  principal  part. 
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consequently,  the  steady  state  equation  is  hyperbolic. 

For  the  special  case  of  a  —  0,  the  rotation  matrix  in  (1.1)  becomes  the  identity,  R  —  /,  and 
(1.1)  has  a  functional  that  is  monotone  decreasing  (subject  to  natural  boundary  conditions), 

E  =  J  |Vu|  dxdy,  ^  =  -  J  [div  (^f|)]2  dx  dy  <  0.  (2.6) 

This  special  case  is  related  to  models  for  geometric  motion  by  mean  curvature  [11],  Gen¬ 
eralizations  of  geometric  evolution  equations  have  been  considered  in  the  field  of  image 
processing  [2,5,8,20].  These  latter  models,  generally  called  Perona-Malik  equations  [24],  are 
based  on  a  different  class  of  generalizations  than  (1.1);  however,  we  will  show  that  they  share 
some  important  properties.  For  (1.1)  with  0  <  a  <  7t/2  there  is  no  decreasing  functional 
corresponding  to  (2.6).  However,  we  shall  identify  a  Liapunov  function  for  one- dimensional 
solutions. 


2.2  The  continuous  one-dimensional  model 


For  any  nonzero  a  —  (oq,  a2)T  6  R2,  the  linear  function 

v(x,  y)  —  aix  +  a2y  (2.7) 

is  a  solution  of  the  steady-state  equation  (2.5),  one  which  describes  uniform  shearing.  In  this 
paper  we  study  only  one-dimensional  perturbations  of  (2.7),  i.e.  solutions  of  the  form 

v(x,  y,  t)  —  a\x  +  a2y  +  w(x,  t).  (2.8) 

Without  loss  of  generality,  throughout  this  paper  we  take  jaj  =  1  in  (2.8),  and  for  simplicity 
we  take  advantage  of  rotational  invariance  of  the  equation  (1.1)  to  make  the  one-dimensional 
perturbation  be  in  the  ^-direction.  Observe  that  for  functions  of  this  form  dyv  =  a2;  thus, 
provided  that  a2  ^  0,  Vu  will  never  vanish,  thereby  avoiding  the  singularity  in  (1.1).  This 
ansatz  is  motivated  partly  by  having  observed  such  one-dimensional  solutions  as  the  large¬ 
time  limit  of  two-dimensional  simulations  of  (1.1)  and  partly  by  our  desire,  in  this  initial 
investigation  of  ill-posed  problems,  to  work  in  a  context  where  much  of  the  behavior  can  be 
derived  rigorously. 

To  examine  ill-posedness  of  the  one-dimensional  problem,  suppose  that  w(x,t )  in  (2.8)  is  an 
exponential,  w  —  e^ix+xt.  Since  Vn0  =  a  and  arg£  =  0,  (2.2)  may  be  simplified  to  show  that 
one-dimensional  perturbations  of  the  steady-state  solution  (2.7)  are  ill-posed  if 

—a  <  arga  <  0.  (2.9) 

For  contrast,  if  —  tt  <  arga  <  —a  or  if  0  <  arga  <  7r  —  a,  then  (2.7)  is  well-posed  in  the 
one-dimensional  context.  In  this  respect,  the  one-dimensional  problem  differs  greatly  from 
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^(^max)  — 

1 

F(S crit)  =  COS  Ct 

/  \ 

1  \ 

- -  F  cos  a 

/ 

y 

F  —  cos  a 

-5  0  5 


s 


Fig.  2.  The  nonlinear  flux  function  F(s)  for  0  <  4>  <  a  (Case  1). 

the  two-dimensional  problem:  (1.1)  is  always  ill-posed  when  plane  waves  in  all  directions  are 
allowed  as  perturbations  of  (2.7). 

On  substitution  of  (2.8)  into  (1.1),  we  obtain 

<™> 

where 


F(s) 


a  +  s  §i  \ 
a  +  s  ei  |  /  ’ 


(2.11) 


here  RT  is  the  transpose  (and  inverse)  of  (2.1),  and  §i  =  (1,0)T.  Recalling  that  \a\  =  1,  we 
define  an  angle  <fi  =  —  arga  (note  the  minus  sign),  so 

a  =  (cos  (j),  —  sin  ^)T;  (2-12) 


without  loss  of  generality  we  restrict  (f>  to  the  range  0  <  <f>  <  tt.  We  shall  consider  (2.10)  on 
the  interval  0  <  x  <  1  and  impose  homogeneous  Dirichlet  boundary  conditions  on  w. 


w(0,t)  =  0,  w(l,t)  =  0.  (2-13) 

In  terms  of  the  two-dimensional  PDE  (1.1),  we  are  seeking  a  solution  on  the  vertical  strip  Q  = 
{{x,y)  e  [0,1]  x  R}  of  the  particular  form  (2.8)  satisfying  boundary  conditions  v(0,y,t)  = 
a^y-,  and  v(l,y,  t)  =  oi  +  a2V ■  It  is  interesting  to  note  that  in  the  cases  (f>  =  0  and  (j)  =  a, 
the  boundary  of  is  characteristic  for  the  steady-state  problem  (2.5).  Incidentally,  the  case 
of  periodic  boundary  conditions  for  (2.10),  i.e.  w(  1  +  x,t)  =  w(x,t),  is  also  covered  by  our 
analysis  apart  from  a  possible  nonzero  mean  value  for  w(x). 

The  properties  of  PDE  (2.10)  depend  on  the  form  of  the  nonlinear  flux  function  F(s).  Using 
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(2.14) 


(2.12),  in  terms  of  cf),  (2.11)  can  be  rewritten 

cos(cx  —  (f>)  +  scoscx 
Vl  +  2s  cos  (j>  +  s2 

Figure  2  shows  a  graph  of  F(s)  for  one  choice  of  a  and  <j>.  Note  that,  for  every  choice  of  cx 
and  <p.  F(s )  is  a  bounded  function  for  all  s,  taking  values  in  the  range  —  cos  a;  <  F(s)  <  1, 
with  the  lower  bound  approached  as  s  — »  — oo.  As  s  — )■  oo, 

.  sin  cx  sin  <5  _9.  , 

F(s)  =  cos  cx  H - H  0(s  2).  (2.15) 

s 

The  upper  limit,  F(s)  =  1,  is  achieved  at  smax,  the  unique  maximum  of  F.  where  the  two 
vectors  in  the  inner  product  (2.11)  are  parallel:  i.e.,  arg(a  +  smaxei)  =  —a.  This  critical  point 
is  given  by 


smax  =  —  sin(cx  —  <^>) /  sin  ex.  (2.16) 

This  point  separates  the  two  intervals  (— oo,smax)  and  (smax,oo),  where  F(s)  is  monotone 
increasing  and  decreasing,  respectively. 


Differentiating  the  one-dimensional  equation  (2.10)  with  respect  to  x.  we  obtain  a  nonlinear 
diffusion  equation  for  the  slope  s(x,t )  =  wx(x,t), 


ds 

dt 


or 


ds 

dt 


d_ 

dx 


(2.17) 


where  D(s )  =  F'(s )  is  the  nonlinear  diffusion  coefficient.  For  s  <  smax,  D(s )  >  0  and 
hence  (2.17)  is  a  linearly  well-posed  nonlinear  forward  diffusion  equation,  while  for  s  >  smax, 
D(s)  <  0  and  (2.17)  is  a  linearly  ill-posed  backward  diffusion  equation  (see  Figure  3b).  We 
observe  that  for  F(s)  given  by  (2.14),  the  corresponding  diffusion  coefficient  is 


D(s)  = 


sin  a  sin  <j> 

(1  +  2s  cos  (j)  +  s2)3^2 


(  S  Smax  )  • 


(2.18) 


The  one-dimensional  Perona-Malik  equation  [20],  wt  =  (p(wl)wx)x,  is  of  the  form  (2.17), 
with  F{wx )  —  p(wx)wx.  For  this  equation,  F'(s)  —  p{s2)  +  2 p'(s'2)s'2  also  becomes  negative 
at  large  s  for  appropriate  functions  p(s2)  of  interest  [20]. 


Equation  (2.10)  has  the  Liapunov  function 

C  =  l  v(wx)dx  where  V(s)  =  J  F(s)ds.  (2.19) 

Subject  to  Dirichlet  or  other  appropriate  boundary  conditions,  the  Liapunov  function  is 
monotone  decreasing,  with  its  evolution  given  by 

^  =  -  Jo  [dxF(wx)f  dx  <  0.  (2.20) 
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Fig.  3.  The  nonlinear  flux  function  F(s)  and  the  corresponding  nonlinear  diffusion  coefficient 
D(s )  =  F'(s)  for  cf>  >  2a  (Case  3). 

Up  to  an  additive  constant,  the  anti-derivative  of  F(s)  in  (2.14)  is  given  by 


We  define  scrit  as  the  value  of  the  slope  in  (— oo,  smax)  where  F(s)  equals  its  limit  as  s  — >  oo: 
be.,  .F (scrit)  =  cos  ct,  or  arg(a  +  scrjtei)  =  —2a:  (see  Figure  2).  This  value  is  given  by 

sCrit  =  —  sin(2o!  —  cf))/sm(2a).  (2.22) 

We  will  show  that  the  long-time  behavior  of  the  discrete  model  of  (2.10)  is  related  to  whether 
Smax  and  scrit  are  positive  or  negative.  We  distinguish  three  cases  as  follows,  indicating  the 
corresponding  relations  between  a  and  < j i,  for  each: 


Case  1  : 

•Scrit  ^  •Smax  ^  0 

0  <  <f>  <  a 

Case  2  : 

Scrit  ^  0  <C  Smax 

cl  <  (j)  <  2a 

(2.23) 

Case  3  : 

0  <C  Scrjt  “s  Smax 

2a  <  <fi  <  7r 

Note  that  in  Case  1,  we  have  F'(0)  <  0  (see  Fig.  2),  so  the  trivial  equilibrium  solution 
w  =  0  is  ill-posed  to  one-dimensional  perturbations.  More  generally,  we  claim  that  any 
possible  smooth  solution  must  be  ill-posed  in  Case  1.  To  see  this,  observe  that  the  boundary 
conditions  (2.13)  force  solutions  to  maintain  zero  average  slope,  f  wxdx  =  0.  Hence  in  Case  1, 
every  non-trivial  continuous  w(x )  will  be  ill-posed  in  some  subset  of  0  <  x  <  1,  where  its 
slope  is  positive.  This  argument  for  the  one-dimensional  problem  (2.10)  independently  re¬ 
derives  and  generalizes  the  earlier  linear  stability  result  (2.9).  In  contrast,  for  Cases  2  and  3 
the  trivial  solution  of  (2.10)  is  well-posed.  Moreover,  for  initial  data  with  max(ryx)  <  smax, 
(2.10)  is  strictly  forward  parabolic.  Consequently  (i)  a  maximum  principle  for  wx  can  be 
established,  (ii)  the  L-2  norm  of  wx  evolves  according  to 

—  w2xdx  =  -  D{wx)w2xxdx<  0,  for  wx  <  smax,  (2.24) 

dt  Jo  Jo 
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and  (iii)  this  class  of  solutions  converge  to  w  —  0  as  t  — >  oo.  While  Case  1  can  be  distinguished 
from  Cases  2  and  3  by  the  respective  ill-  or  well-posedness  of  w  —  0,  we  will  show  that  Cases 
2  and  3  differ  in  whether  or  not  w  —  0  is  the  unique  equilibrium  solution  of  the  discretized 
version  of  (2.10). 


2.3  The  semi-discrete  one- dimensional  problem 


As  discussed  in  the  introduction,  we  consider  a  continuous-time/discrete-space,  approxima¬ 
tion  to  (2.10,  2.13).  Specifically,  given  a  positive  integer  N  >  2,  let  Ax  —  1/N  be  the 
uniform  spacing  of  grid  points  in  a  finite-difference  discrete  solution  wn(t)  «  w(nAx,t),  for 
n  —  0, 1,  2, . . .  ,  N.  The  discrete  solution  evolves  according  to  the  coupled  system  of  ( N  —  1) 
nonlinear  ordinary  differential  equations  for  n  —  1,  2,  ...N  —  1, 


dwn  _  1  f  /wn+1-wn\  /wn-wn_1\) 

dt  ~  Ad  l  Ax  J  V  Ax  J J  ’ 

where  the  boundary  conditions  (2.13)  become 

w0  —  0,  wN  —  0. 


(2.25) 


(2.26) 


For  N  —  2,  (2.25,2.26)  reduces  to  a  single  first  order  equation  for  w p  similarly  for  N  —  3, 
the  general  model  reduces  to  an  autonomous  phase  plane  system  for  W\  and  w2.  Recently 
and  independently,  work  on  these  low-dimensional  models  was  done  in  [35]  for  the  study  of 
clustering  instabilities  in  granular  gases.  Our  analysis  of  the  solutions  of  (2.25,  2.26)  applies 
for  all  N  >  2,  but  our  focus  will  be  the  consideration  of  large  values  of  N  and  the  continuum 
limit,  N  — »  oo. 

For  brevity,  we  will  write  the  arguments  of  F(-)  in  (2.25)  as 


w 


71+1/2  — 


Wn+i  ~  Wr 


Ax 


(2.27) 


This  is  a  second-order-accurate  centered  finite-difference  approximation  of  the  spatial  deriva¬ 
tive,  i.e.  w'n+1/2(t)  =  wx((n+^)Ax,  t)+0( Ax2).  This  discrete  system  has  a  Liapunov  function 
analogous  to  (2.19), 


JV-l 


£(™n)  =  V(W'n+ 1/2)  AX. 

71=0 


(2.28) 


Using  summation  by  parts,  it  can  be  shown  similarly  that  the  discrete  Liapunov  function  is 
monotone  decreasing, 
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Formally,  the  discrete  model  (2.25)  converges  to  the  PDE  (2.10)  as  Ax  —$■  0  up  to  0( Ax2) 
errors.  Retaining  terms  through  order  0(Ax4),  we  derive  the  effective  PDE  for  (2.25), 


dw 

dt 


d_ 

dx 


F(wx)  + 


Ax2  1  d 
24  wxx  dx 


+  0{  Ax4). 


(2.30) 


Linearizing  this  equation  about  w  =  0,  we  obtain 


dw 

dt 


m 


( d2w  Ax2  d4uA 

+  12  dx4  ) 


+  0{  Ax4). 


(2.31) 


For  D( 0)  <  0,  equation  (2.31)  contains  a  destabilizing  second-order  term  and  a  regularizing 
fourth-order  term.  This  balance  of  terms  occurs  in  many  models  of  physical  systems,  includ¬ 
ing  the  Cahn-Hilliard  equation  for  binary  mixtures  [7,29]  and  the  Kuramoto-Sivashinsky 
equation  [15]  from  combustion  theory,  image  processing  models  [20,8,2],  biological  systems 
[23,21],  and  spatially  discrete  mechanical  systems  [28,26].  Solutions  of  these  equations  exhibit 
regions  where  the  solution  is  smooth,  which  are  separated  by  thin  layers  with  large  gradi¬ 
ents;  for  large  times  these  layers  eventually  merge  together  [29,30].  We  observe  analogous 
coarsening  behavior  in  solutions  of  (2.25). 

For  -D(O)  >  0  the  fourth  order  linearized  modified  PDE  (2.31)  is  high-frequency  unstable. 
However,  as  noted  in  the  introduction,  the  discrete  finite  difference  model  (2.25)  eliminates 
the  dynamics  of  all  length  scales  smaller  than  the  fundamental  grid-spacing.  Consequently 
the  spatially  discrete  system  (2.25)  is  not  subject  to  such  unbounded  growth  rates  at  high 
frequencies. 

Since  F(s)  is  bounded,  solutions  of  (2.25)  exist  for  all  time  and  grow  at  most  linearly  in 
time.  We  show  that,  in  fact,  independent  of  Ax  =  1/7V,  the  solution  is  bounded  for  all  time 
in  the  max  norm:  i.e.  there  exists  a  C  such  that 


max  I wn(t)  I  <  C. 

n  1  v 

In  view  of  the  boundary  condition  (2.26),  boundedness  of  wn(t)  is  an  immediate  consequence 
of  the  following  proposition. 

Proposition  2  If  w(t)  is  a  solution  of  (2.25)  and  if  Sq  =  min„  w'n+1/2( 0),  then  for  all  n  and 
all  t  >  0 


w'n+ 1/2(*)  >  min(scrit,50). 


(2.32) 
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Note  that  if  w(x,t )  were  a  twice  differentiable  solution  of  (2.10),  then  w  would  satisfy  the 
stronger  estimate  wx(x,t)  >  min(smax,  S0)  where  S0  =  infs  wx(x,  0).  Indeed,  since  F(s)  is 
monotone  increasing  on  (— oo,smax),  the  maximum  principle  gives  this  estimate.  However, 
because  of  the  discretization,  solutions  of  (2.25)  do  not  satisfy  this  stronger  estimate.  Nev¬ 
ertheless,  the  derivation  of  (2.32)  is  modeled  on  the  proof  of  the  maximum  principle. 

Proof:  For  brevity  let  s*  denote  the  RHS  of  (2.32).  Obviously,  (2.32)  is  initially  satisfied  at 
time  t  =  0.  Suppose  that  for  all  grid  points  and  for  all  t  <  t*  we  have  w'(t)  >  s*,  and  suppose 
that  for  some  grid  point  n,  we  have  w'n+1/2 (t*)  =  s*.  Then,  taking  differences  of  (2.25),  we 
calculate  that  at  this  point  and  t  =  t*, 

=  {F(W'n+ 3/2)  -  2F(s*)  +  F(w'n_  1/2))  . 

We  claim  that 

F(s*)  <  F(w'n+ 3/2) 

and  similarly  for  Ff ir’(J_1/2),  from  which  it  follows  that  at  f  =  t*. 

dwn+l/2  > 

dt  t=u  ~ 

To  prove  the  claim,  observe  that  s*  <  w'n+^i2.  If  w'n+3j2  <  scrjt ,  then  (2.33)  follows  from 
the  fact  that  F(s)  is  monotone  on  an  interval  containing  (— oo,scrit).  On  the  other  hand,  if 
Sent  <  <+3/2 >  then 

FM  <  ^(Scrit)  <  ^(<+3/2) » 
the  latter  inequality  in  (2.35)  being  apparent  from  Figure  2. 

By  itself,  condition  (2.34)  is  not  sufficient  to  exclude  the  possibility  that  w'n+1  ^2(t)  <  s*  for 
some  t  >  t*.  However,  as  in  the  proof  of  the  maximum  principle  [18],  we  may  derive  (2.32) 
from  the  above  argument  by  taking  the  limit  of  slightly  modified  functions  for  which  (2.34) 
becomes  a  strict  inequality.  ■ 


3  The  discrete  steady-state  solutions 

The  trivial  solution,  wn  =  0,  is  an  equilibrium  of  (2.25).  In  this  section  we  determine  the 
nontrivial  equilibrium  solutions  of  (2.25).  We  show  that  the  number  of  solutions  depends  on 
(f),  corresponding  to  the  three  cases  identified  in  (2.23). 


(2.35) 


(2.33) 


dw'n+ 1/2 
dt  t=u 
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3.1  Equilibrium  solutions  of  the  discrete  model 


If  an  (N  +  l)-component  vector  {wn}  is  an  equilibrium  solution  of  (2.25,  2.26)  then  there  is 
a  constant  F  such  that  the  slopes  (2.27)  all  satisfy 

F(w'n+ 1/2)  =  F,  n  =  0,l,.. .  ,N-1.  (3.1) 

From  the  boundary  conditions  (2.26),  the  N  values  of  the  discrete  slope,  {w'n+ 1/2\,  are  also 
subject  to  the  global  constraint 


E  <+i/2=  0-  (3-2) 

n= 0 

This  condition  is  the  discrete  analogue  of  f  wxdx  =  0,  the  consequence  of  the  homogeneous 
Dirichlet  boundary  conditions  (2.13).  Thus,  we  have  identified  equilibrium  solutions  {wn}  of 
(2.25,  2.26)  with  solutions  ({w'n+1/2},  F)  of  (3.1,  3.2). 

Considering  (3.1)  alone,  we  note  that  nontrivial  solutions  are  possible  only  when  F(s )  =  F 
has  multiple  solutions.  If  cos  a  <  F  <  1,  there  are  two  values  of  the  slope  s,  call  them  si 
and  s2,  such  that  F(s)  =  F, 


F(s  0  =  F(s2)  =  F.  (3.3) 

From  the  properties  of  F(s)  given  by  (2.14)  it  is  clear  that  Si  and  s2  must  bracket  smax. 
When  F  «  1  both  roots  are  close  to  smax;  as  F  decreases  the  roots  separate  continuously. 
As  F  —f  cos  a;,  one  solution  of  (3.3)  becomes  large,  Si  00,  while  the  other  one  approaches 
the  limit  s2  —t  scrjt.  Using  (2.14)  we  can  eliminate  F  in  (3.3)  to  express  Si  in  terms  of  s2, 

sin(2q— 2<j>)  _ 

Sl  =  H(s2)  =  - —  (3.4) 

+>crit  S2 

This  formula  shows  that,  as  F  varies,  (si,  s2)  traces  out  a  portion  of  a  hyperbola,  as  shown 
in  Figure  4.  This  hyperbola  is  necessarily  invariant  under  the  interchange  of  Si  and  s2.  Its 
horizontal  and  vertical  asymptotes  are  sij2  — >■  scrit-  Changing  the  values  of  a  and  <j>  affects  the 
position  of  this  hyperbola  in  the  plane  and  changes  the  number  of  solutions  of  the  discrete 
problem  (3.1,  3.2). 

To  complete  the  description  of  the  nontrivial  equilibrium  solutions,  we  now  examine  the 
constraint  (3.2).  To  enumerate  the  equilibrium  solutions,  we  define  I\  as  the  number  of  grid 
positions  where  the  slope  u4+1/2  is  sq:  the  remaining  N  —  K  positions  will  have  slope  s2. 
Consequently  (3.2)  reduces  to 


Ks  1  +  (N  -  K)s2  =  0. 


(3.5) 
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Si 


Fig.  4.  The  construction  for  the  weak  and  strong  (s®  ,  s®)  solutions  in  Case  2,  given  by 

the  intersection  points  of  the  hyperbola  (3.4)  and  the  line  (3.5). 


Case  3’ 

Case  3’ 

/ 

Case 

_ _ -^'^Ca se  2’ 

Case  1’ 

Fig.  5.  (left)  The  upper  bound  for  Case  2',  <f>  <  /3(a,  N.  1),  for  N  =  10,  20, 40,  •  •  •  ,  and  the  continuum 
limit,  (3  — »  2a  as  N  — >  oo,  and  (right)  /3(a,N,K)  for  K  in  1  <  K  <  N  with  fixed  N. 

Graphically,  finding  the  intersection  points  of  the  hyperbola  (3.4)  and  the  line  (3.5)  in  the 
(si,s2)  plane  yields  all  of  the  nontrivial  equilibrium  solutions  at  given  (q;,^),  see  Figure  4. 
Eliminating  sq  from  (3.4,  3.5)  yields  a  quadratic  equation  for  s2.  This  equation  has  two  real 
solutions  if  the  resulting  discriminant  is  positive, 

(N  —  2 K)2  sin2(2o;  —  (j))  +  AK(N  —  K )  sin(2o;)  sin(2o;  —  2(f))  >  0.  (3.6) 


The  case  of  equality  in  (3.6)  corresponds  to  the  case  when  the  line  (3.5)  is  tangent  to  the 
hyperbola  (3.4)  and  defines  an  upper  bound  for  (ft  in  Case  2.  We  write  this  bound  as  the 
function  [3,  (f>  =  /3(a1  N,  K)1  in  terms  of  a,  N ,  K  with  a  <  (3  <  2a  (see  Figure  5). 

Using  /3(a1N1I\)1  we  can  define  three  cases  for  the  equilibrium  solutions  of  the  discrete 
problem  (2.25)  that  correspond  to  the  cases  given  by  (2.23).  In  the  limit  that  N  — >  oo,  for 
any  fixed  K1  we  find  that  the  upper  bound  for  Case  2  is  <f>  =  (3  — >  2a1  corresponding  to 
(2.23)  (see  Figure  5).  For  finite  N,  the  condition  (ft  =  f3(a,  TV,  K)  defines  the  degenerate  case 
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Fig.  6.  The  bifurcation  diagram  for  the  weak,  strong,  and  trivial  equilibrium  solutions. 

where  the  quadratic  for  s2  has  a  double  root.  For  every  (j>  <  (3  there  are  two  distinct  real 
roots  for  s2,  which  we  will  call  the  strong  and  weak  solutions  (denoted  by  superscripts  S  and 
W),  s®  <  see  Figure  6.  For  the  range  0  <  (j>  <  a  corresponding  to  Case  1,  (see  (2.23)),  sf 
and  s ^  have  opposite  signs,  see  Figure  6.  At  (j)  —  a ,  the  weak  solution  becomes  degenerate, 
s ^  =  0,  and  it  intersects  the  branch  of  trivial  solutions,  w  =  0.  As  we  will  describe  later, 
this  intersection  point  is  a  transverse  bifurcation  that  is  connected  with  a  change  in  stability 
of  these  solution  branches.  For  the  range,  a  <  (j>  <  (3,  corresponding  to  Case  2  (as  N  — »  oc) 
both  strong  and  weak  solutions  yield  negative  slopes  for  s2.  Finally,  for  the  range  (f  >  A 
corresponding  to  Case  3  (as  N  — >  oc),  there  is  only  the  trivial  solution  w  =  0. 

We  summarize  these  results  on  the  classification  of  nontrivial  equilibrium  solutions  in  terms 
of  s2  as: 

For  any  N  >  2  and  for  each  K  with  1  <  I\  <  N: 

Case  V  :  Two  solutions,  s|<0<s^  0  <  <^>  <  ck 

Case  2'  :  Two  solutions,  sf<s^<0  a  <  (f>  <  / 3(a ,  N,  K)  (3.7) 

Case  3'  :  No  nontrivial  solutions  ft  (a,  N ,  K)  <  (j>  <  n. 

Case  V  is  the  same  as  Case  1  of  (2.23);  Cases  2'  and  3'  approach  Cases  2  and  3  of  (2.23) 
as  N  — >  oc,  for  which  j3  — v  2a.  The  corresponding  si  values  are  given  by  (3.5),  si  — 
—  (N  —  K)s2/K.  For  given  values  of  N  and  K ,  there  are  distinct  equilibrium  solutions 
{ wn }  resulting  from  different  spatial  arrangements  of  the  two  slopes.  Clearly  s\  and  s2  must 
have  opposite  signs,  and  if  K  <C  N  then  sfj  |s2|,  for  example  see  Figure  7.  These  solutions 
with  K  N  can  be  described  as  having  jump  discontinuities  at  the  grid  points  with  the 
larger  slope  si,  with  the  magnitude  of  the  jump  being  siAx  =  Si/N:  the  remainder  of  the 
solution  is  piecewise  linear  with  slope  s2.  We  will  refer  to  si  as  the  “jump”  or  “shock”  slope, 
and  to  s2  as  the  “background”  or  “mean-held”  slope,  respectively. 

At  this  point  we  note  that  the  above  description  actually  provides  a  redundant  enumeration 
of  the  set  of  equilibrium  solutions.  This  is  due  to  double  counting  of  the  solutions  caused  by 
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Fig.  7.  Two  equilibria  for  Case  V  with  a  =  7t/8,  <fi  =  a/2  and  N  =  100;  (a)  a  strong  shock  solution 
with  K  =  2  jumps,  and  (b)  a  much  smaller-amplitude  weak  shock  solution  for  K  =  1. 

the  symmetry  under  the  interchange, 


si  S2  and  K  N  —  K,  (3.8) 

i.e.  if  {si,  S2,  K,  N  —  k}  describes  a  solution  of  (3.3)  and  (3.5),  then  so  does  {s2,  si,  N  —  K,  K}. 
To  eliminate  this  redundancy,  let  us  dehne  the  two  slopes  in  a  nontrivial  equilibrium  solution 
with  one  slope,  s+,  that  is  greater  than  smax,  s+  >  smax,  and  the  other  less  than  smax, 
s_  <  smax.  For  enumeration  of  the  solutions,  we  can  associate  s+  with  si  and  henceforth 
define  K  as  the  number  of  grid  cells  with  slope  s+, 

Ks+  +  (N-  K)s _  =  0  (3.9) 

This  description  of  the  two  equilibrium  slopes  will  be  of  particular  importance  for  the  later 
analysis  of  the  stability  and  dynamics  of  solutions  of  (2.25).  This  is  because  s+  >  smax 
corresponds  to  the  ill-posed  regime  for  the  PDE  (2.17)  with  D(s)  <  0,  while  s_  <  smax 
corresponds  to  well-posed  behavior  with  D(s )  >  0.  Except  in  the  special  case  of  the  Case 
V  weak  solution,  the  descriptions  of  solutions  based  on  the  ill-posed  and  well-posed  slopes, 
(s+,s_),  are  equivalent  to  the  descriptions  in  terms  of  the  jump  and  background  slopes, 
(si,S2). 

We  now  derive  asymptotic  estimates  for  the  two  families  of  weak  and  strong  solutions  with 
a  single  jump,  K  =  1,  in  the  limit  of  large  N .  In  this  limit,  s+  tends  to  infinity,  and  from 
Figure  2  we  see  that  s_  tends  to  scrit;  more  precisely,  by  condition  (3.5)  we  obtain 

SS  =  Sent  +  OiN-1)  <  0,  (3.10) 

4  =  -(iV-lKrit+O(l)>0. 

We  will  call  these  solutions  with  slopes  (s®  ,  s® )  strong  shocks.  They  correspond  to  the  lower 
branch  of  solutions  shown  in  the  bifurcation  diagram  Figure  6,  where  we  note  that  (j>  =  0 
yields  scr;t  =  —  1,  i.e.  the  endpoint  of  that  branch  of  the  bifurcation  diagram. 

For  K  =  1  as  N  — »  oo,  the  other  solution  is  given  by 
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(3.11) 


sw  =  -s0/(N-1)+O(N-2)1 

sw=  S0  +  O(lV -1), 

where  s o  is  the  nonzero  slope  such  that  F(sq)  =  F( 0),  specihcally, 

So  =  —  sin(2o;  —  2  (f))/ sin(2o;  —  ).  (3-12) 

We  will  call  these  solutions  weak  shocks.  As  N  — )■  oo,  these  {wn}  solutions  are  small  in 
amplitude,  they  scale  as  0(N~1)  -»  0  everywhere,  with  the  size  of  the  jump  being  s^Ax  ~ 
sq/N  (see  Figure  7).  The  weak  shocks  correspond  to  the  upper  branch  of  solutions  in  Figure  6. 
At  4>  =  a,  so  =  0,  and  as  noted  earlier,  the  weak  shock  solution  coincides  with  the  trivial 
solution  wn  =  0. 

Note  that  similar  results  for  multiple-jump  weak  and  strong  shock  solutions  can  be  derived 
for  any  fixed  I\  =  1,2,3,...  in  the  limit  that  N  — >-  oo.  Then,  given  the  values  of  the  two 
equilibrium  slopes,  s+  and  s_,  and  their  spatial  distribution,  say  the  values  n^,  with  k  = 
1,2,  ...,/F,  of  the  grid  points  where  the  slopes  s+  =  w'r jfc+1/2  occur,  then  the  solution  {wn} 
can  be  reconstructed  explicitly  from 

n 

wn  =  J2wj+ 1/2  Ax,  n  =  0,  1,2*....,  N.  (3.13) 

l=o 


3.2  Comparison  of  discrete  equilibria  with  generalized  solutions  of  the  PDE  (2.10) 


As  described  above,  the  equilibrium  solutions  of  the  discrete  problem  (2.25)  are  piecewise 
linear  functions  with  K  finite  jumps.  For  the  case  of  K  =  1  jump,  applying  the  bound¬ 
ary  conditions  (2.26)  and  summing  the  difference  quotient  wn+ 1/2,  (2.27),  over  n  yields  the 
equilibrium  solution 


wn  =  s-  [nAx  -  H{n  -  [n*  +  |])]  ,  (3-14) 

where  Ax  =  1  /TV,  Ft(-)  is  the  Heaviside  function,  and  n*  is  the  value  of  n  for  which  w'n+1/.2  = 
s+.  For  N  — >-  00,  this  solution  is  the  discrete  analogue  of  a  weak  solution  of  the  PDE  (2.10). 

Formally,  an  equilibrium  solution  of  (2.10)  has  dxF(wx )  =  0,  or  equivalently  F(wx )  =  F. 
If  F  =  cos  a  then  there  is  a  family  of  piecewise  linear  weak  equilibrium  solutions  with  a 
single  finite-jump  discontinuity.  The  mean-held  equilibrium  slope  is  the  finite  solution  of 
F{wx )  =  cosct,  that  is  wx  =  scrit,  (2.22).  Consequently  we  can  write  the  equation  for  an 
equilibrium  solution  with  a  single  discontinuity  as  wx  =  scrit[l  —  S(x  —  2*)],  where  5  is  the 
Dirac  delta  function  and  0  <  x*  <  1  is  the  shock  position.  Integrating  this  and  applying 
(2.13)  yields  the  steady-state  solution 

w(x)  =  sciit[x  -  H(x  -  x^)].  (3.15) 


17 


This  one-parameter  family  of  solutions,  parametrized  by  the  shock  position  in  the  domain 
0  <  x*  <  1,  is  the  continuum  limit  of  (3.14),  since  as  N  — >  oo,  s+  — >  oo  and  s_  — >  scrjt. 

However,  (among  many  other  solutions)  piecewise  linear  weak  equilibrium  solutions  can  be 
constructed  with  any  countable  number  of  positive  jump  discontinuities.  We  are  content  to 
mention  in  passing  the  formal  similarity  between  solutions  of  the  discrete  system  and  gener¬ 
alized  solutions  of  the  ill-posed  PDE  (2.10).  More  definitive  statements  about  the  solutions 
of  the  PDE  require  careful  analysis  [14,20].  In  the  following  sections,  we  study  in  detail  the 
stability,  local  instabilities  and  global  dynamics  of  the  discretized  model  (2.25). 


4  Stability  of  the  equilibria 


The  following  result  on  the  stability  and  long  time  evolution  of  solutions  of  (2.25)  for  the 
cases  defined  in  (3.7)  will  be  proved  in  this  section. 

Proposition  3  On  the  stability  of  equilibrium  solutions  of  (2.25): 

Case  V  Only  strong  shock  solutions  of  (2.25)  with  a  single  jump  discontinuity  are  stable. 
Case  2!  Strong  shock  solutions  with  a  single  discontinuity  and  the  trivial  zero  solution  are 
the  only  stable  equilibria. 

Case  3'  The  zero  solution  is  stable  and  in  fact  globally  attracting. 

In  Cases  V  and  2'  for  almost  every  initial  condition  (i.e.,  not  on  the  stable  manifolds  of  the 
unstable  equilibria),  every  solution  approaches  a  stable  equilibrium  as  t  — oo. 


4-1  The  Liapunov  function:  a  necessary  condition  for  stability 


We  now  make  use  of  the  Liapunov  function  (2.28) 


=  E  V(wn+ 1/2)  Axi 
n= 0 

to  establish  some  fundamental  stability  results  for  the  discrete  system  (2.25).  Recall  from 
(2.29)  that  £  is  monotone  decreasing  as  w  evolves.  Indeed,  £  defines  (2.25)  as  a  gradient 
system  in  the  form 


dwn  1  0£ 

dt  Ax  dwn  ’ 


n  =  1,2, •••IV-  1. 


(4,1) 


In  particular,  w  is  an  equilibrium  solution  of  (2.25)  if  and  only  if  it  is  a  critical  point  of  £. 
£  is  bounded  from  below  and  tends  to  infinity  as  \w\  — >  00.  Therefore  £  must  always  have 
at  least  one  local  minimum  corresponding  to  a  linearly  stable  solution. 
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In  Case  3',  w  =  0  is  the  only  critical  point  of  C ;  thus  from  the  properties  of  C,  w  =  0  must  be 
a  minimum  of  C.  Therefore,  in  Case  3',  the  trivial  solution  is  stable  and  globally  attracting. 


Turning  to  Cases  1'  and  2',  the  following  result  eliminates  most  candidates  in  the  search  for 
stable  equilibria. 

Proposition  4  If  w°  e  RN+1  is  an  equilibrium  solution  of  (2.25)  derived  from  (3.3,  3.5) 
with  K  >2,  then  w°  is  not  a  local  minimum  of  C. 


Proof:  Consider  an  equilibrium  solution  w°  with  s+,  s_  given  by  (3. 4, 3. 5)  for  a  given  K  >  2. 
To  demonstrate  that  such  a  solution  is  not  stabl,  we  show  that  C  is  not  a  minimum  at  w° 
by  differentiating  along  an  appropriate  curve  of  vectors  in  R  v+1,  {w (<?)},  through  w°.  We 
construct  this  one-parameter  family  of  near-equilibrium  states  by  perturbing  two  of  the  K 
values  where  the  slope  is  w'n+1/2  =  s+  (at  points  n  =  rq  and  n  =  n2)  then  the  finite  difference 
quotients  w'n+1/2(q )  satisfy 


s+  +  q 

if 

n  —  ni, 

s+  —  q 

if 

n  =  n2, 

S+ 

if 

n  —  nk  for  k 

S- 

otherwise, 

(4.2) 


where  we  note  that  for  q  —  0,  we  recover  the  equilibrium,  w( 0)  =  w°,  while  the  constraint 
(3.2)  is  satisfied  for  all  q.  Then  the  Liapunov  function  is 

£(»(«))  =  4  [r(S+  +  q)  +  V(s+  -q)  +  ( K  -  2)V(s+)  +  (N  -  K)V(s-)l .  (4.3) 

At  q  —  0  the  first  derivative  of  C(w(q))  vanishes,  and  the  second  derivative  satisfies 

=  vv/"(s+)  =  vF<*+)  <  °- 

where  the  final  inequality  follows  from  s+  >  smax-  Consequently  we  conclude  that  (4.2)  with 
q  —  0  is  not  a  local  minimum  of  C  and  hence  w°  with  K  >  2  is  an  unstable  equilibrium.  ■ 


d2c 

dq2 


4-2  Linear  stability  analysis 


Having  used  the  Liapunov  function  to  establish  the  instability  of  all  equilibria  with  more 
than  one  jump,  we  turn  to  linear  stability  to  analyze  the  trivial  solution  and  single  jump 
(K  —  1)  equilibria. 
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First,  at  the  trivial  solution  wn  =  0,  the  linearization  of  (2.25)  is 

dw  T  171/  /r\\  wn-\-i  -  2u;n  +  wn-i 

~di  =  Lw^F  (0) - - ' 

where  L  =  N2F'(0)T  is  the  (N  —  1)  x  (TV  —  1)  symmetric  tridiagonal  matrix  with 

-2  1 
1  -2  1 
1  -2  1 


1 


The  matrix  T  is  the  standard  finite  difference  centered  second  derivative  operator  with  ho¬ 
mogeneous  Dirichlet  boundary  conditions  [16].  Since  T  is  negative  definite,  with  eigenvalues 
Oj  =  —  4sin2(i7rj/Ar)  for  j  =  1,  2,  ..N  —  1,  the  trivial  solution  of  (2.25)  is  stable  if  F'{ 0)  >  0 
—  in  Cases  2'  and  3'  —  and  is  unstable  if  F'{ 0)  <0  —  in  Case  V . 

To  complete  the  stability  analysis,  we  examine  the  linearization  of  (2.25)  at  single  jump 
(I\  =  1)  weak  and  strong  shock  equilibria  for  Cases  T,  2'.  We  will  summarize  the  results 
below  in  Proposition  6.  Note  that  from  Proposition  4,  we  can  already  eliminate  the  possibility 
that  the  single-jump  weak  shock  solution  (s)^,  s^),  is  a  stable  equilibrium  in  Case  1',  since 
from  (3.8)  it  maps  onto  a  solution  with  s+  >  smax  with  K  =  N  —  1  >  2,  however  we  will 
mention  it  in  the  discussion  for  completeness. 

The  linearization  of  (2.25)  may  be  analyzed  in  terms  of  matrices  containing  two  blocks 
similar-in-form  to  (4.5).  To  facilitate  the  calculation,  we  introduce  the  notation  for  the 

matrix  (4.5);  the  subscript  of  course  specifies  the  dimension,  and  the  superscripts  refer  to  the 
Dirichlet  boundary  conditions  at  both  end  points  of  the  interval.  Extending  this  notation, 
we  shall  write  T ^  for  the  analogous  M-dimensional  operator  with  Dirichlet  boundary 
conditions  at  the  left  endpoint  and  Neumann  boundary  conditions  at  the  right  endpoint: 
i.  e..  the  M  x  M  matrix  with  rows  as  in  (4.5)  except  with  the  last  row  replaced  by 

(0  ...  0  1  -1).  (4.6) 

Similarly,  is  obtained  by  modifying  the  first  row  of  T^td\  and  by  modifying  both 
the  first  and  last  rows. 

For  an  equilibrium  with  one  jump  located  at  grid  point  ni  =  /,  that  is  ^/+1/2  =  s+;  the 
linearization  of  (2.25)  may  be  viewed  as  a  perturbation  of  a  block  diagonal  matrix, 

T„  =  diag(T<*),TK_1).  (4.7) 


(4.4) 
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Specifically,  the  linearization  is  given  by  the  symmetric  tridiagonal  operator 


L(e)  =  A2F'(s_)(T0-eP), 


where 

r_  F'(s+) 
F'(s-Y 


and  the  perturbation  P  is  given  by 

P  =  diag(Oj_i,  A,  0 jv-/-2), 

with  0 m  denoting  the  M  x  M  zero  matrix  and  A  the  2x2  matrix 


— i  A 
i-i  / 


(4.8) 


(4.9) 


(4.10) 


(4.11) 


Note  that  P  has  three  blocks,  the  2x2  middle  block  overlapping  the  corners  of  the  two 
blocks  in  (4.7).  Here  we  have  assumed  that  2  <  /  <  N  —  3;  the  cases  with  a  jump  adjacent 
to  either  endpoint,  which  are  simpler,  are  left  for  the  reader. 


Proposition  5  T0  —  eP  is  negative  definite  if  —  oo  <  e  <  (N  —  1)  1  and  has  one  positive 
eigenvalue  if  e>  (TV  —  l)-1. 


Proof:  We  shall  show  that 


det[— (T0  -  eP)]  =  1  -  (N  -  l)e. 

As  shown  by  Givens  (see  [16]),  the  determinant  dN- 1  of  an  (TV  —  1)  x  (N  —  1)  symmetric, 
tridiagonal  matrix  with  entries  M  =  {my}  is  generated  recursively  by  dj  =  mjjdj- 1  — 
(■ m,j-ij)2dj-2  for  j  =  2,  3, ...,  N  —  1  with  d0  =  1  and  d\  =  mu-  In  applying  this  algorithm  to 
M  =  —  (T0  —  eP),  for  j  =  2, . . .  ,  I  —  1  and  for  j  =  I  +  2, . . .  ,  N  —  1  the  recursion  relation  is 

dj  —  2dj-i  dj-2, 

while  for  the  two  values  of  j  between  these  two  ranges 

di  =  (1  —  e)di- 1  —  di- 2  and  dI+ 1  =  (1  —  e)di  —  e2d/_i. 

In  the  first  range  of  j,  from  the  recursion  relation,  we  obtain  dj  =  1  +  j  >0  for  j  = 
2,3,  •••,/—  1.  In  the  second  range  of  j,  after  applying  the  special  cases  for  j  =  /  and 
j  =  1+1  given  above,  we  find  dj  =  1  —  ej  for  j  =  /,  •  ■  •  ,  N  —  1.  Thus,  the  determinant  is 
given  by  dN- 1  =  1  —  (N  —  l)e,  as  claimed. 

Further,  from  Givens’  theorem  [16],  the  number  of  positive  eigenvalues  of  T0  —  eP  is  given 
by  the  number  of  sign  changes  in  the  sequence  {dj}.  Note  that  T0  is  negative  definite,  hence 
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for  e  =  0  there  are  no  positive  eigenvalues.  Since  {dj}  for  j  >  I  is  monotone  decreasing, 
only  a  single  sign  change  can  occur,  if  dN- 1  <  0.  Since  dN- i(e)  =  0  has  a  simple  zero  for 
e  =  1/(N  —  1),  the  matrix  T0  —  eP  has  a  single  positive  eigenvalue  for  e  >  1/(7V  —  1).  ■ 

Proposition  6  The  linear  stability  of  the  single  jump  (K  =  1 )  equilibria  breaks  down  into 
cases  as  given  by  (3.7): 

Case  V  :  The  strong  shock  solution  is  stable,  while  the  weak  shock  is  highly  unstable  with 
N  —  1  positive  eigenvalues. 

Case  21  :  The  strong  shock  is  stable,  while  the  weak  shock  solution  is  unstable  with  one 
positive  eigenvalue. 

Proof:  We  begin  with  the  strong  shock,  (s®,s®).  It  can  be  shown  that  ss_  <  smax  for 
0  <  (j>  <  f3(a,  N,  1)  and  for  any  N.  Therefore  the  factor  F'(s-)  in  (4.8)  is  always  positive  for 
strong  shocks.  To  estimate  e  =  ef  in  (4.9),  we  use  the  asymptotic  form  (3.10)  for  (s®  ,  s+)  in 
the  formula  (2.18)  for  D(s)  =  F'(s )  : 


(1  +  2scrit  cos  +  sgrit)3/2 

(®max  ^crit )  ^crit 


=  0(N~ 2)  >  0. 


(4.12) 


Consequently  since  e  =  0(N  2)  <C  (N  —  1)  1  as  N  — >-  oo,  by  Proposition  5,  all  of  the 
eigenvalues  of  L(e)  are  negative  and  the  strong  shock  is  stable  for  Cases  V  and  2'. 

In  contrast,  for  the  weak  shock  (sw,  s^),  in  the  limit  N  — >  oo,  we  find  from  (2.18)  and  (3.11) 
that 


e  =  ef  ~ - g°  5max - ^  =  0(1)  >  0.  (4.13) 

Smax(l  +  2scrit  COS  (j)  +  S2rit)'  1 

Therefore,  for  the  weak  shock,  e  =  0(1)  >>  (N  —  l)-1  as  N  — »  oo,  and  T0  — eP  is  not  negative 
dehnite  but  has  one  positive  eigenvalue  for  both  Cases  V  and  2! .  For  Case  2r,  sff  <  smax,  so 
F'(sff)  is  positive  and  the  weak  shock  is  unstable  with  one  positive  eigenvalues.  For  Case  T, 
sff  >  Smax,  so  the  multiplicative  factor  F'(sff)  is  negative  and  the  IV  —  1  negative  eigenvalues 
of  T0  —  eP  become  N  —  1  unstable  positive  eigenvalues  for  L(e)  for  the  weak  shock.  ■ 

As  was  shown  above  in  Proposition  6,  for  Case  V,  since  it  is  stable,  the  strong  shock  must 
correspond  to  a  minimum  of  the  Liapunov  function  C.  From  (2.28),  this  value  of  C  is  inde¬ 
pendent  of  the  position  of  the  jump  within  the  domain.  Therefore  all  N  of  the  single-jump 
strong  shock  solutions  are  stable.  The  same  argument  can  be  applied  for  the  strong  shock  in 
Case  2' .  These  results  are  independent  of  the  linearized  analysis  done  in  this  section.  Linear 
stability  analysis  shows  that  the  position  of  the  shock  and  the  influence  of  the  boundary  con¬ 
ditions  do  weakly  effect  the  values  of  the  eigenvalues,  but  they  do  not  change  the  stability 
of  the  solutions. 
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Fig.  8.  Convergence  of  different  numerical  methods  for  an  initial  value  problem  in  Case  V  for  system 

(2.25) . 

5  Representative  dynamic  simulations 

In  this  section  we  present  the  results  of  a  representative  set  of  numerical  simulations  of  the 
dynamics  of  system  (2.25).  These  simulations  illustrate  some  of  the  differences  in  behavior 
in  the  three  cases  in  Proposition  3.  The  simulations  also  guide  the  analysis  of  the  nonlinear 
dynamics  in  the  following  sections. 

We  begin  with  a  brief  discussion  of  the  numerical  methods  used  for  the  simulations.  As  was 
discussed  above,  while  the  continuum  PDE  (2.10)  is  ill-posed,  for  any  finite  N  the  discrete  sys¬ 
tem  (2.25)  is  well  posed,  with  global  solutions.  Subject  to  typical  analytic  constraints  [16,17], 

(2.25)  can  be  solved  numerically  using  any  appropriate  method  for  integrating  systems  of 
coupled  ODEs.  Figure  8  shows  the  results  of  a  convergence  study  as  the  discrete  time-step 
approaches  zero,  At  — >  0,  for  a  typical  initial  value  problem.  We  tested  several  standard  ex¬ 
plicit  and  implicit  schemes.  For  sufficiently  small  At  all  of  the  methods  showed  convergence 
with  the  expected  order  of  accuracy  (see  Figure  8).  The  implicit  midpoint  method,  written 
in  general  form  as 

v.m+l  _  ...m 

"  At  ~  =^.(sK‘+1+M,m>)’  (5.1) 

where  w™  «  wn(tm )  and  tm  =  mAt,  had  the  smallest  error  coefficient  and  was  used  for  all 
of  the  following  numerical  simulations. 

In  Figure  9  we  consider  the  evolution  of  the  solution  for  an  initial  value  problem  in  Case  1', 
with  a  —  7t/8,  <j)  —  a/2,  and  wn( 0)  =  sm(nn/N)  for  n  —  0,1,2,---  ,  JV  with  N  —  100. 
As  described  above,  in  Case  V ,  the  stable  steady  state  is  piecewise  linear  with  a  single 
jump.  The  intermediate  dynamics  leading  up  to  this  state  are  rather  complicated.  Figure  9a 
shows  the  initial  unstable  behavior;  the  solution  rapidly  develops  a  large  number  of  jump 
discontinuities,  forming  what  is  sometimes  called  a  “staircase  pattern.”  This  stage  of  the 
evolution  can  be  compared  to  spinodal  decomposition  of  binary  mixtures  [7],  where  large 
numbers  of  phase  interfaces  develop  from  an  unstable  initial  state.  Figure  9b-d  shows  the 
generic  mode  of  evolution  for  longer  times;  the  sizes  of  the  jump  discontinuities  evolve  on 
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Fig.  9.  Evolution  for  a  Case  1'  initial  value  problem;  (a)  short  time  evolution  up  to  time  ta  cap¬ 
tures  the  unstable  phase,  where  smooth  initial  data  forms  lots  of  discontinuities,  (b-d)  coarsening 
behavior  at  times  tf,,tc  leading  to  the  final  steady-state  with  a  single  shock,  see  t,/. 

slower  time  scales.  This  regime  will  be  described  as  coarsening  dynamics,  where  most  of  the 
phase  interfaces  collapse  leaving  larger  intervals  where  the  mean  field  holds. 

Due  to  the  global  constraint  (3.2),  while  some  of  the  jumps  grow,  others  must  decay.  This 
behavior  is  illustrated  in  a  different  form  in  Figure  10b,  where  the  values  for  all  of  the 
slopes,  w'n+1/2(i),  n  —  0 ...TV  —  1  are  plotted  as  functions  of  time.  From  this  figure  we  note 
that  at  any  time  t  0.1  there  are  a  small  number  of  points  with  large  positive  slopes 
(jumps  with  s+  >  smax),  while  most  of  the  grid  points  have  a  small  negative  slope  (the 
mean  field  with  s_  <  smax)-  From  Figure  10a  we  observe  that  the  dynamics  for  (2.25)  has  a 
monotone  decreasing  Liapunov  function  (2.28)  that  experiences  a  sequence  of  rapid  declines 
coinciding  with  the  collapse  of  each  successive  jump.  Ultimately,  the  single-jump,  stable, 
strong  shock  is  approached  for  <  t  — ? >  oo  (see  Figure  9d).  Incidentally,  the  location  of  the 
final  jump  depends  very  sensitively  on  the  initial  data  and  on  the  simulation  parameters. 
In  a  series  of  numerical  experiments  like  that  of  Figure  9,  but  with  perturbed  initial  data 
wn( 0)  =  (1+e)  sin(7 xn/N)  with  e  =  O(10-12),  we  found  that  even  such  tiny  perturbations  lead 
to  discontinuous  changes  in  the  position  of  the  steady-state  jump.  This  extreme  sensitivity 
is  a  clear  reflection  of  the  ill-posed  nature  of  the  underlying  problem. 

For  contrast  with  Figure  9,  we  present  the  evolution  of  the  problem  in  Case  3',  with  (f)  —  2.1a , 
starting  from  the  same  initial  conditions,  see  Figure  11.  Figure  11a  shows  that  for  long  times, 
the  solution  converges  to  the  stable  trivial  solution,  w  —  0.  However,  since  the  initial  data  is 
partially  ill-posed,  with  the  slopes  of  the  initial  condition  satisfying  w'n+1j2  >  smax  for  some 
range  in  x,  a  staircase  pattern  composed  of  finite  jumps  with  large  slopes  develops  for  t  & 
0.1.  For  comparison  with  Figure  10b,  Figure  lib  shows  the  more  complicated  intermediate 
dynamics  for  the  large  slopes  in  Case  3',  before  they  all  decay  to  zero. 
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Fig.  10.  (a)  The  evolution  of  (a)  the  Liapunov  function,  and  (b)  the  values  of  the  local  slopes  w'n+1j,2 
for  n  =  0, 1,  2,  ...N  —  1  plotted  as  a  function  of  time  for  the  Case  V  problem  in  Fig  9. 


0  0.5  1  0.01  0.1  1  10 

xn  t 


Fig.  11.  Evolution  for  the  initial  value  problem  in  Case  3;;  (a)  wn  profiles  and  (b)  slopes  w'n+1/2  as 
a  function  of  time. 

As  described  in  Proposition  3,  for  Case  2'  both  the  trivial  solution  w  —  0  and  the  K  =  1 
strong  shock  solutions  are  stable.  One  way  to  illustrate  this  bi-stability  is  to  plot  the  Liapunov 
function  C  for  the  one-parameter  family  of  piecewise  linear  functions  with  I\  —  1  (see  Fig.  12), 

w'ni+l/2  =  S+i  Wn+ 1/2  =  for  71  ±  nl-  (5‘2) 

Then,  as  a  function  of  the  slope  s+  at  the  jump,  the  Liapunov  function  (2.28)  takes  the  form 

jC(s+)  =  T  [v'(s+)  +  (Af  -  1)V  (-^Tt)]  -  (5'3) 

This  is  a  double  well  potential  with  minima  corresponding  to  the  trivial  state  w  —  0  and  the 
strong  shock  solution  (see  Fig.  12a).  The  weak  shock  is  an  unstable  equilibrium  corresponding 
to  a  maximum  of  C  and  separates  the  basins  of  attraction  of  the  two  stable  states  (see 
Fig.  12b).  While  this  description  is  quite  suggestive,  it  is  also  somewhat  misleading  for  the 
dynamic  evolution.  This  is  because  solutions  of  (2.25)  starting  from  initial  data  given  by 
(5.2)  do  not  remain  within  the  family  (5.2)  with  s+  =  s+(t).  Nevertheless,  this  discussion 
serves  to  point  out  the  significance  of  the  unstable  weak  shock  as  a  boundary  for  the  basins 
of  attraction  of  the  trivial  and  strong  shock  solutions.  The  existence  of  bi-stability  in  Case  21 
in  a  model  of  clustering  of  granular  gases  [35]  of  the  same  form  as  (2.25)  with  N  —  3  was 
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Fig.  12.  (a)  The  double-well  Liapunov  function  and  (b)  equilibria  for  single-jump  (K  =  1)  piecewise 
linear  solutions  in  Case  2'. 


Fig.  13.  The  basin  of  attraction  for  the  uniform  state  w  =  0  for  a  <  4>\  <  </>2  <  4>s  <  ■  ■  •  <  /3 (ok,  IV,  1) 
in  Case  2'  given  in  terms  of  the  two-parameter  family  of  initial  conditions  (5.4). 


studied  in  connection  with  hysteretic  effects. 

Next,  we  illustrate  the  dependence  on  the  initial  conditions  of  the  large-time  limit  of  the 
solution  in  Case  2'.  Consider  the  discrete  problem  (2.25)  with  N  =  100  and  a  two-parameter 
family  of  small-amplitude  initial  data  (ei,  e 2  small), 

w(x,  0)  =  ei  sin(7rx)  +  e2  sin(27nr).  (5.4) 

In  Case  2',  the  basin  of  attraction  of  the  trivial  solution  depends  on  the  value  of  (j>,  with 
a  <  (j)  <  ft  (a,  N,  1).  We  plot  the  boundary  of  the  basin  of  attraction  of  the  trivial  solution 
w  —  0  for  various  values  of  (j)  in  this  range,  see  Figure  13.  All  initial  conditions  within  the 
basin  converge  to  w  —  0,  large-amplitude  solutions  outside  the  boundary  converge  to  the 
stable  strong  shock  solution.  In  the  limit  (j)  — >  a,  the  basin  shrinks  to  the  origin,  as  in  Case 
1,  where  w  —  0  is  unstable  and  the  shear  band  is  the  global  attractor.  In  the  limit  <j)  — >  9, 
the  basin  of  attraction  for  the  strong  shock  solution  shrinks  to  a  point,  where  there  is  a 
saddle-node  bifurcation  and  the  weak  and  strong  shocks  coalesce  (see  Fig.  6).  The  trivial 
solution  w  —  0  then  remains  as  the  unique  (global)  attractor. 
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Fig.  14.  Illustration  of  the  1/3  scaling  laws  during  the  coarsening  regime  of  evolution;  (a)  K/N , 
the  density  of  jumps  in  the  solution,  s+  >  smax,  and  (b)  the  slope  of  the  average  jump  .s/v's  for  a 
set  of  simulations  with  a  range  of  values  for  N  =  100.. 25600. 

6  Scaling  laws  and  considerations  of  the  continuum  limit 


In  this  section  we  study  in  detail  the  coarsening  dynamics  of  solutions  of  (2.25)  in  Case  V 
with  initial  conditions 


w(x,  0)  =  w(x)  =  Hsin(7nc), 


(6.1) 


with 

0  <  \A\  <  Hmax, 

where  Hmax  =  |smax|/7T-  For  this  range  of  A.  the  condition 

dxw(x)  >  smax,  (6.2) 

is  satisfied  everywhere  and  coarsening  dynamics  is  observed  throughout  the  entire  interval, 
0  <  x  <  1.  Recall  from  Section  2.2  that  in  Case  1,  the  trivial  solution  w  =  0  is  linearly 
ill-posed  and  small-amplitude  initial  data  satisfying  (6.2)  will  be  unstable  everywhere  in 
0  <  x  <  1.  Grid-scale  instabilities  develop  everywhere  and  evolve  to  a  staircase  profile  (see 
Fig.  9ab)  followed  by  coarsening  dynamics,  during  which  the  jumps  vanish  one-by-one  until 
just  one  remains.  It  is  easier  to  analyze  the  dynamics  for  small  initial  data  than  the  simulation 
shown  in  Fig.  9,  where  A  =  1.  We  will  consider  initial  value  problems  with  small  initial  data 
for  Case  V  of  (2.25)  to  avoid  complications  that  may  be  introduced  by  bi-stability  in  Case  2' 
and  non-uniform  spatial  instabilities  for  large  amplitude  data.  While  we  show  results  for  the 
specific  example  with  initial  data  with  A  =  10-9  in  Case  V  with  a  =  7t/8  and  <p  =  7t/16, 
and  smax  ~  —0.855,  for  simulations  with  N  =  2m  x  100  for  m  =  0,1,2,...,  8  grid  points, 
simulations  with  other  data  strongly  suggest  that  the  features  of  the  ensuing  evolution  are 
universal,  for  sufficiently  large  N,  for  all  initial  data  satisfying  (6.2). 

In  section  3,  we  defined  equilibrium  jumps  as  grid  points  where  the  discrete  slope  was  larger 
than  smax,  w'n+1/2  >  smax.  Here,  we  apply  the  same  criterion  for  counting  the  number  of 
jumps,  K(t),  in  an  evolving  solution  {wn(t)}.  The  remaining  N  —  K  grid  points,  where 
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w'n+1j2  <  «max,  are  called  the  background.  Examination  of  the  simulations  shows  that,  for 
large  times,  K(t),  the  number  of  finite  jumps  at  time  t ,  satisfies  the  scaling  law 

/  t \-V3 

m  ~  C  (-J  ,  (6.3) 

for  some  constant  C.  A  remarkable  collapse  (even  for  short  times)  of  the  data  from  simula¬ 
tions  with  different  values  of  N  occurs  if  we  re-express  (6.3)  as  a  scaling  law  for  K/N ,  the 
density  of  jumps  in  the  interval, 


~  C(NH)-113.  (6.4) 

This  scaling  behavior  is  exhibited  in  Figure  14a,  which  contains  data  from  simulations  with 
N  =  2m  x  100  for  m  =  0, 1,  2, ...,  8  grid  points. 

A  direct  consequence  of  (6.4)  is  a  scaling  law  for  the  average  slope  for  a  jump  at  time  t. 
To  see  this  consequence,  we  note  that  the  global  constraint  on  the  slopes  (3.2)  holds  for  all 
times.  Let  s+g  and  savg  be  the  averages  at  time  t  of  the  jump  (s  >  smax)  and  background 
( s  <  smax)  slopes  respectively.  Then  in  terms  of  these  averages,  (3.2)  yields  a  generalization 
of  (3.9)  valid  for  all  times, 


Ks+g  +  (N-  K)saAg  =  0.  (6.5) 

respectively.  For  long  times,  savg  «  scrit  and  K  decreases,  leading  to  K  <C  N,  so  (6.5)  reduces 
to 

4>s~-^SCT„~-^i(A'2t)1/3.  (6.6) 

Figure  18b  shows  that  for  large  times,  the  average  slope  at  a  jump  does  follow  this  scaling 
behavior. 


A  useful  heuristic  picture  of  the  solution  may  be  extracted  from  formula  (6.3),  considered 
at  a  fixed  time  t0,  as  N  — >■  oo.  In  the  coarsening  dynamics  regime,  the  solution  {wn(t0)}  is 
approximately  piecewise  linear  in  n,  with  intervals  of  width  0(K~l)  =  O^to/N)1^3)  where 
the  slope  is  close  to  scrit  alternating  with  single  grid-cells  with  large  negative  slopes  of  order 
0((N2to )2/3).  As  suggested  by  Figure  9a,  wherever  staircase  pattern  forms,  it  evolves  so  that 
{wn(t0)}  continues  to  follow  the  initial  data  w(x)  in  some  approximate  or  locally-averaged 
sense.  Based  on  this  observation  we  introduce  two  norms  to  investigate  how  the  solution 
evolves  from  the  initial  data. 


Specifically,  we  consider  two  L2  norms,  for  the  depature  of  wn(t)  from  the  initial  condition 
w(x),  (6.1),  and  for  the  departure  of  the  slopes  w'n+1j2(t)  from  w^ix), 


Ew(t,N ) 


F  \wn(t)  —  w(nAx)\2  Ax 


1/2 


(6.7) 
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Fig.  15.  Evolution  of  the  norms,  Ew  and  Es  for  a  series  of  simulations  in  the  limit  N  — >  oo;  (a)  very 
short-time  behavior  showing  linear  growth  and  IV-dependent  linear  instability,  and  (b)  evolution 
on  the  fast  timescale,  showing  saturation  of  the  phase  separation  instabilities. 

1  1/2 

\w'n+i/2(t)  ~  w^{nAx)  I2  Ax  .  (6.8) 

.  n 

Figure  15  shows  plots  of  these  norms  for  a  series  of  simulations  with  resolutions  N  —  2m  x  100 
points  with  m  —  0, 1,  2,  •  •  •  ,8.  From  Figure  15a,  we  see  that  for  very  short  times  Ew(t)  shows 
slow  exponential  growth  with  rate  Ai  ~  —tt2F'(0).  This  is  to  be  expected  since  (6.1)  is  a 
multiple  of  the  lowest  order  eigenvector  of  the  linearization  (4.5).  However,  since  this  problem 
is  ill-posed,  this  smooth  evolution  is  very  quickly  overwhelmed  by  strongly  unstable  high 
frequency  modes  that  generate  grid  oscillations.  These  instabilities  also  grow  exponentially, 
with  rates  on  the  order  of  the  largest  eigenvalue,  A^v-i  =  0(N2)  as  N  -*  oo  (see  Figure  15a). 
This  instability  can  be  regarded  as  the  initial  stage  of  phase  separation  -  the  formation  of 
large  gradients  in  the  solution.  Due  to  the  maximum  principle,  Proposition  2,  these  grid 
oscillations  cannot  grow  indefinitely,  but  saturate  and  lead  to  another  stage  of  dynamics. 


Es(t,N)  = 


From  the  growth  of  the  strong  instabilities,  we  note  that  the  timescale  of  the  dynamics  at 
the  end  of  the  regime  of  linearized  growth  is  t  —  0(N~2).  In  fact,  this  timescale  holds  for 
all  of  the  longer  time  dynamics  of  the  solution.  Figure  15b  shows  that  for  longer  times,  both 
norms,  Ew  and  Es  depend  on  the  rescaled  time,  r  =  N2t.  In  fact,  there  is  a  remarkable 
collapse  of  the  results  from  all  of  the  simulations  onto  limiting  curves  independent  of  N  that 
hold  after  the  instabilities  have  saturated.  Figure  15b  shows  that  for  large  r,  the  norms  very 
closely  follow  the  power-law  scaling  with  exponent  1/3, 

Ew(t,N)  ~  ^(!V2t)1/3,  E2(t,N)  ~  Cs(N2t)1/3.  (6.9) 

Changing  the  point  of  view,  fixing  t  and  letting  N  — »■  oo,  from  the  scaling  of  Es(t),  we 
note  that  the  maximum  slope  in  the  solution  will  diverge  like  0(N 2/3)  (see  also  (6.6)),  see 
Figure  16a.  Further,  note  that  the  scalings  for  norms  (6.9)  at  a  fixed  time  simplify  to 

Ew(t,  N )  =  0(N -1/3)  E2(t,  N )  =  0(iV2/3).  (6.10) 


From  this  we  observe: 
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Fig.  16.  Scaling  laws  for  the  properties  of  the  solution  at  a  fixed  time,  t  =  1,  in  the  limit  that 
N  — »  oo;  (a)  the  evolution  norms  =  0(iV2/3)  and  Ew  =  C^iV-1/3),  (b)  K/N ',  the  density  of 
jumps  in  the  solution. 

•  For  any  fixed  positive  time  t,  as  N  — >  oo,  in  terms  of  the  L 2  norm,  the  solution  wn(t)  will 
not  have  evolved  from  the  initial  data,  since  Ew  =  0(Ar_1/3)  — >  0. 

•  For  any  fixed  positive  time  t,  as  N  — »•  oo,  in  terms  of  the  H 1  norm,  the  solution  blows  up, 
since  E 2  =  0(N 2</3)  — >  oo. 

These  two  observations  for  the  finite-time  behavior  of  the  continuum  limit,  iV  ->  oo,  describe 
a  solution  that  evolves  from  the  initial  data  only  by  instantaneously  developing  infinitesimally 
small  jump  discontinuities.  This  singular  behavior  for  the  continuum  limit  shows  that  the 
discreteness  of  model  (2.25)  is  essential  to  its  dynamics  for  any  finite  N. 

Note  that  re-writing  (6.9)  for  Ew  in  the  form 

/M1/3 

Ew(t,N)~Cw(- J  ,  (6.11) 

shows  that  (2.25)  has  a  very  long  timescale,  t  =  0(N),  associated  with  the  coarsening 
dynamics.  This  is  the  timescale  that  describes  the  very  slow  overall  evolution  of  the  system. 
As  described  earlier,  during  coarsening,  the  large  number  of  jumps  initially  created  during 
the  initial  phase  separation  regime  will  systematically  collapse  producing  successive  solutions 
with  fewer,  larger- amplitude  jumps.  In  the  final  section,  we  present  some  analysis  of  this 
dynamic  behavior. 


7  Intermediate  Dynamics:  Coarsening 

The  dynamic  simulations  of  the  previous  sections  suggest  that  while  we  have  thoroughly 
studied  the  steady  states  and  asymptotic  stability,  a  complete  understanding  of  the  behavior 
of  (2.25)  requires  an  examination  of  the  complicated  intermediate  dynamics  of  the  system 
as  well.  For  ill-posed  initial  data,  i.e.  for  data  with  max(vnx)  >  smax,  the  formation  of  a  large 
number  of  jumps  in  the  solution  creates  very  unstable  intermediate  states.  From  Proposi- 
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tion  4,  we  know  that  there  are  no  stable  states  with  more  than  K  —  1  jumps.  Consequently, 
the  dominant  feature  of  the  evolution  for  all  finite  times  will  be  a  type  of  coarsening  dy¬ 
namics  —  a  process  continually  reducing  the  number  of  jumps  in  the  solution  until  a  stable 
steady  state  is  achieved.  In  this  section  we  will  use  two  approaches  to  examine  the  dynamic 
behavior  at  a  single-step  transition,  from  K  to  K  —  1  jumps.  First,  we  use  linear  analysis 
to  study  the  unstable  equilibria,  then  we  also  consider  an  approximate  reduction  of  the  full 
system  to  a  lower  dimensional  nonlinear  system. 


7.1  Linearization  at  two-jump  equilibria 


To  assess  the  transient  timescale  for  the  collapse  of  a  jump  discontinuity,  we  estimate  the 
positive  (unstable)  eigenvalue  of  the  linearization  of  (2.25)  at  an  equilibrium  with  two  strong 
shocks.  As  was  done  in  Section  4,  the  linearization  can  be  analyzed  as  a  perturbation  of  a 
matrix  with  block  structure,  but  now  with  three  blocks.  Specifically,  if  the  jumps,  separated 
by  L  grid  points,  are  at  the  grid  points  rq  =  I  and  n2  =  I  +  L,  then 

L(e)  =  Ar2F,(s_)(T0  —  eP), 

where 

T0  =  diag(Tjdn),  T^n\  tJIm), 

e  is  given  by  (4.9),  and 

P  =  diag(0/_i,  A,  0L_2,  A,  On-i-l-2)-  (7.3) 

Here  A  is  the  2x2  matrix  given  by  (4.11)  In  the  following  calculation  we  shall  let  N  — »  oo 
while  keeping  I/N  and  L/N  fixed. 

The  matrix  T0  in  (7.2)  is  negative  semidefinite  with  one  zero  eigenvalue  associated  with  the 
eigenvector 

u  =  (0/,  eL,  0n-i-l-i)T ,  (7.4) 

where  0/  denotes  the  /-dimensional  zero  vector  and  €  RL  is  given  by  =  ( 1,  1, . . .  ,  1)T. 
This  eigenvector  spans  the  kernel  of  the  second-difference  operator  with  Neumann  boundary 
conditions,  T^.  To  leading  order  in  perturbation  theory, 

Amax  -  —eN2F'(s-)  (7.5) 

\u,u) 

where  (•,  •)  denotes  the  normalized  Euclidean  inner  product  on  RA'_1,  (u,v)  —  J2unvn  Ax. 
It  is  readily  computed  that  (u,  Pw.)  =  —2/N  and  that  (u,  u)  —  L/N ;  thus, 

2  eN2F'(ss_) 

^max  ^  7 


(7.2) 


(7.1) 
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Fig.  17.  Transient  evolution  for  multiple-jump  solutions;  (a)  time-evolution  for  the  collapse  of 
a  jump  in  a  K  =  2  solution  (solid  line)  compared  with  the  linear  growth  rate  (7.9)  (dashed 
line),  (b)  comparison  of  the  maximum  unstable  eigenvalue  for  IT-jump  equilibria  (dots)  with  the 
perturbation-theory  estimate  given  by  (7.11)  (dashed  line). 

To  estimate  e,  we  use  the  analogue  of  (3.10)  for  a  two-jump  (K  =  2)  strong  shock  equilibrium, 


A  max  1 


(7.11)- 


ss_  =  sCTit  +  0(N~1)  <0,  (7.7) 

4  =  — (iV  —  2Krit/2  +  0(l)  >  0. 


Using  (4.9)  and  (7.7),  we  obtain  the  value  of  e  for  the  two-jump  strong  shock  solution  in 
terms  of  the  previous  result  (4.12), 


_  _  _  F'iS+)  /|J3  _  /-u /U-2a 

e  -  e2  -  ^  ~  4ex  -  <J(N  ), 


Consequently,  we  obtain  the  estimate, 2 


(7.8) 


Amax 


8  sin  a  sin  (p 

T  e2  ' 


(7.9) 


Of  course  1/Amax  defines  the  time  scale  for  the  collapse  of  a  two-shock  meta-stable  equi¬ 
librium;  in  particular,  since  L  scales  like  N/K,  the  collapse  time  is  proportional  to  N  as 
iV  — >  oo.  Figure  17a  illustrates  the  accuracy  of  this  linear  estimate  for  the  evolution  for  the 
collapsing  jump,  starting  from  a  small  perturbation  of  the  K  =  2  strong  equilibrium  (7.7), 
see  Figure  7a.  The  linear  growth  rate  (7.9)  gives  a  very  good  estimate  of  the  evolution  until 
the  jump  has  almost  completely  collapsed,  w'ni+1/2  ~  smax. 

More  generally,  this  argument  can  be  extended  to  show  that  if  there  are  K  jumps  at  points 
rii,  n2,  •  •  •  ,  uk  where  I\  <C  N,  then  the  perturbation-theory  estimate  for  the  largest  positive 

2  Note  that,  like  e,  the  smallest  eigenvalues  of  To  are  0(N~2).  As  a  check  on  the  accuracy  of 
applying  perturbation  theory  when  the  perturbations  are  of  the  same  order  as  the  eigenvalues,  we 
showed  that  the  second  order  correction  to  the  eigenvalue  is  0(eiV-2),  smaller  than  the  first  order 
correction  by  a  factor  of  N.  Moreover,  we  determined  the  lowest  eigenvalue  of  T0  —  eP  using  Sturm 
sequences  (see  [16]),  and  this  yielded  the  same  results  as  perturbation  theory. 
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eigenvalue  of  the  linearization  is 


A 


max 


r\,/ 


2  esKN2F'(ss_) 
min k(nk  -  nk- 1)' 


(7.10) 


If  the  jumps  are  approximately  equally  spaced,  then  the  denominator  is  Lm jn  ~  N/K.  For 
solutions  with  I\  jumps,  ~  J\2ef,  where  ef  =  0(N~ 2),  see  (4.12).  Consequently,  our  crude 
estimate  for  the  timescale  for  the  transition  from  K  to  K  —  1  jumps  is 


1  s 2  ■  N 

x  ^cnt  1 V 

-  rsj - # 

Amax  2  sin  a  sin  <f>  K 3 


(7.11) 


A  comparison  of  this  estimate  of  the  maximum  growth  rate  with  the  largest  unstable 
eigenvalue  of  the  i7-jump  strong  equilibria  is  shown  in  Figure  17b,  which  confirms  that 
Amax  =  0(K3).  For  K  >  2,  the  strong  shock  equilibria  have  K  —  1  closely  spaced  positive 
eigenvalues,  and  a  simplified  single-mode  linearized  analysis  may  be  insufficient  to  describe 
the  dynamic  evolution  of  the  problem  for  moderate  to  larger  values  of  K. 


We  note  that  the  results  of  the  linearized  analysis,  (7.11),  suggest  that  the  scaling  law  for 
K(t)  should  have  the  exponent  one  half  rather  than  the  actual  value  of  one  third,  observed 
from  the  simulations  in  Section  6.  The  failure  of  this  linear  estimate  indicates  that  the 
dynamic  solution  wn(t)  does  not  come  arbitrarily  close  to  the  unstable  77-jump  equilibria 
during  the  coarsening  process.  In  consequence,  the  rate  of  collapse  of  the  jumps  is  faster 
than  the  estimate  from  linear  theory. 


7.2  The  jump-diffusion  model 


In  this  subsection,  we  formulate  a  model  that  seeks  to  isolate  the  evolution  of  the  finite 
jumps,  the  main  feature  of  the  nonlinear  coarsening  dynamics.  Let  sn  denote  the  discrete 
slope 


sn  —  wn+ 1/2 


Wn+ 1  ~  Wr, 

Ax 


n=  0, 1,  2,  -  -  -  ,N  -  1. 


(7.12) 


Taking  differences  of  (2.25),  we  write  the  equivalent  slope-evolution  equations, 


dSn  _  F{sn+ 1)  -  2 F(sn)  +  F{sn- 1) 
dt  ~  Ax 2 

Similarly,  at  the  edges  of  the  domain,  n  =  0  and  n  =  N  —  1, 


(7.13) 


dsp  _  F(si)  -  F(fo)  ds jy-i  _  FjsN-i)  -  F(s± v_2) 

dt  Ax2  ’  dt  Ax2 


(7,14) 


Equations  (7.14)  are  derived  from  the  Dirichlet  boundary  conditions  (2.26).  In  particular, 
the  equation  for  sn- i  is  obtained  by  differentiating  the  boundary  condition  constraint  (3.2), 
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in  the  form 


N- 2 

Sn-1  =  -  Sm  (7-15) 

n= 0 

with  respect  to  t  and  using  (7.13)  to  collapse  the  sum. 

Note  that  equations  (7.13),  are  a  second-order  accurate  finite  difference  discretization  of  the 
slope  diffusion  PDE  (2.17): 


1=1^4 

The  jump  in  wn  is  expressed  by 

wn+i  wn  —  snAx.  (7.17) 

This  quantity  corresponds  to  a  jump  in  w  as  Ax  — >  0  only  when  sn  is  large,  sn  =  0(N ) 
as  N  — >  oo.  In  the  intermediate  dynamics,  we  find  that  only  large  positive  values  of  sn  can 
occur,  with  sn  >  smax,  as  was  the  case  for  the  stable  equilibrium  solutions  found  in  section  3. 
Consequently,  as  was  done  in  Section  6,  points  in  the  solution  {wn}  where  the  slope  satisfies 
sn  >  *max  will  be  referred  to  as  jumps,  and  the  remaining  points  with  small  slopes,  sn  <  smax, 
will  be  called  the  background. 

In  Figure  18a,  we  show  a  typical  numerical  simulation  with  N  =  500,  at  a  time  when  there 
are  K  —  20  jumps  on  0  <  x  <  1.  The  evolution  of  a  jump  is  controlled  by  the  second 
difference  52F(sn )  =  F(sn+ 1)  —  2 F(sn)  +F(sn- 1),  which  appears  as  the  numerator  of  (7.14). 
Specifically,  whether  a  jump  will  grow  or  decay  depends  on  if  52F(sn )  is  positive  or  negative, 
respectively.  Figure  18a  shows  a  portion  of  the  solution  wn,  0  <  n  <  180,  containing  seven 
jumps  (label  them  ai,  cr2,  ■  ,a7)  and  Figure  18b  shows  the  corresponding  flux  F(sn).  Note 
that  the  flux  is  continuous  everywhere,  and  the  first  difference,  SF(sn )  =  F(sn+ 1)  —  F(sn ), 
is  piecewise  constant  with  jumps  corresponding  to  those  of  wn.  Moreover,  Figure  18b  shows 
that  the  second  difference  S2F  is  negative  for  jumps  a2  and  <r6,  hence  those  jumps  will  decay 
in  amplitude.  Also  note  that  the  jump  labelled  a4  is  close  to  equilibrium,  since  locally  F(sn ) 
is  nearly  linear  and  hence  S2F(a4)  is  near  zero  (see  Figure  18b). 

Consider  a  solution  of  (7.13,  7.14)  starting  with  K  interior  jumps  at  some  initial  time. 
Numerical  simulations  suggest  that  the  background,  with  sn  <  smax,  equilibrates  on  a  fast 
time-scale,  and  thereafter  evolves  quasi-statically,  being  driven  by  the  slower  evolution  of 
the  jumps.  In  particular,  since  the  background  is  at  quasi-static  equilibrium,  we  deduce  that 
F(sn )  must  be  approximately  linear  between  jumps,  as  indicated  in  Figure  18b.  We  further 
observe  in  this  figure  that  F(sn)  appears  continuous  at  the  jumps.  Consequently,  the  flux  of 
solution,  F(sn ),  is  determined  everywhere  by  the  values  of  F(sn)  at  the  jumps  alone.  With 
these  observations,  we  can  formulate  a  closed  system  of  equations  for  the  evolution  of  the 
jumps  that  is  decoupled  from  the  quasi-static  evolution  of  the  background. 
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Let  ak  =  snk  denote  the  slope  at  the  kth  jump.  Taking  F(sn)  to  be  linear  in  n  between  jumps 
yields  the  model 

F(sn)  =  F(ak)  +  F^Gk+l^ — (n  _  nfe)  for  nk  <  n  <  nk+1.  (7.18) 

Tlk  ■  1  k 

To  specify  this  model  for  the  flux  to  the  left  of  the  first  jump,  and  to  the  right  of  the  Kth 
jump,  define  <r0  =  So  and  uk+ i  =  sjv-i-  Then  F(sn)  is  defined  for  all  0  <  n  <  N  —  1  by 
(7.18)  in  terms  of  ak  for  0  <  k  <  I\.  Since  the  jumps  are  assumed  to  be  interior,  both  s0 
and  sjv-i  are  less  than  smax,  are  part  of  the  background,  and  hence  evolve  quasi-statically. 
In  particular,  the  rates  of  evolution  for  a0  and  ok+i  in  (7.14)  are  of  lower  order  than  those 
of  the  rates  of  evolution  for  the  jumps  {ak}.  Balancing  terms  in  (7.14)  forces  the  conditions 
that  F(s i)  =  F(sq)  and  F(sm~-2 )  =  F(sn- i)-  Consequently,  the  slope  of  F(sn)  adjacent  to 
each  boundary  is  zero,  and  hence 

F(a0)  =  F(ai),  F(aK+1)  =  F(aK) 

After  some  manipulation,  substitution  of  (7.18)  into  (7.20)  leads 

dok  _  1  (F{ok+1)  -  F(ok)  F(ak)  -  F(ak_1)\ 

dt  Ax2  \  nk+ 1  -  nk  nk  -  nk-X  )  ’ 

and  similarly  including  (7.19)  gives  equations  for  o\  and  ok  '■ 

doi  _  F(a2)  -  F(a1)  daK  _  F(aK)  -  . 
dt  (n2  —  ni) Ax2  ’  dt  (uk  —  n r 

This  system  is  similar  in  form  to  (7.14),  but  it  represents  a  vast  reduction  of  the  problem 
when  K  <C  N  —  we  have  to  consider  only  K  coupled  equations  at  the  jumps,  rather  than 
N  equations  at  all  of  the  points  in  the  domain.  Whereas  (7.13)  describes  a  finite  difference 
scheme  for  (7.16)  on  a  uniform  grid,  (7.20)  is  a  discretization  of  (7.16)  on  a  non-uniform 
grid,  given  by  the  positions  of  the  jumps,  {nk}. 

It  is  perhaps  worth  considering  the  evolution  of  jumps  in  the  context  of  the  nonlinear  diffusion 
PDE  (7.16)  for  the  slope  field  s(x ,  t ).  If  the  solution  w(x,  t)  of  (2.10)  contains  a  finite  number 
of  jump  discontinuities  at  locations  xk ,  then  the  slope,  s  =  dxw,  is  a  distribution  containing 
delta- functions  at  xk.  Assuming  the  jump  locations  xk  are  stationary  (independent  of  t). 
we  find  that  dts  also  has  delta  functions  at  the  jumps.  Consequently,  interpreting  (7.16)  in 
the  sense  of  distributions,  we  find  that  the  flux  F(s)  is  continuous  in  space  and  dxF(s )  is 
piecewise  continuous  in  space  with  jumps  at  xk.  Moreover,  at  the  jumps,  we  have 

dt[w]k  =  [dxF(s)]kl 

where  [w]k  =  w(xk  , t)  —  w(xk,t )  denotes  the  kth  jump  in  w.  and  the  right  hand  side  denotes 
the  corresponding  jump  in  dxF(s).  Thus,  in  the  continuum  version  of  the  jump-diffusion 
model  (7.20),  it  is  not  clear  how  to  express  [w]k  in  terms  of  s.  so  the  system  is  not  closed. 


(7.19) 

to  the  equations 
2<k<K-  1,  (7.20) 


F{oK- i) 


r_i  )Ax2 


(7.21) 
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Fig.  18.  (left)  A  portion  of  a  solution  {wn(fo)}  with  N  =  500  containing  K  =  20  jumps,  seven  of 
which  are  on  0  <  x  <  0.35.  (right)  A  plot  of  the  corresponding  values  of  the  flux  function,  F(sn) 
indicating  the  decaying  jumps,  where  52F(s)  <  0. 

This  is  in  contrast  with  the  discrete  model,  for  which  jumps  in  wn  are  related  to  sn  through 
equation  (7.17). 

We  now  justify  the  claim  of  separation  of  time-scales  in  the  discrete  model  when  K  <C  N. 
As  described  earlier,  if  the  slopes  {sn}  for  the  background  correspond  to  a  smooth  function 
s(x,t )  <  smax,  then  as  N  — »  oc,  then  (7.13)  converges  to  the  nonlinear  forward-diffusion 
equation  (7.16)  with  time-scales  independent  of  N.  Hence  t  =  0(1)  for  the  evolution  of 
transients  in  the  background.  In  contrast,  for  large  N,  the  jumps  are  07 .  =  0(N),  and  by 
using  the  asymptotics  of  the  flux  for  large  s  -*  00,  (2.15),  in  (7.20)  yields  the  slow  time-scale 
for  the  evolution  of  the  jumps  as  t  =  0(N).  This  observation  was  obtained  in  section  6.1  for 
the  decay  rate  of  unstable  jumps  in  the  near-equilibrium  case,  (7.11). 

We  now  briefly  review  the  details  of  what  happens  as  a  jump  collapses.  Decay  of  a  jump  occurs 
on  the  0(N )  long  time-scale  if  $2F(ok)  <  0;  the  slope  decreases  to  a  value  with  a *  <  smax. 
Local  equilibrium  will  be  achieved  when  52F{pk)  =  0.  As  the  evolution  proceeds,  the  diffusion 
coefficient,  D{pk )  =  F'(ok )  changes  sign  from  negative  to  positive,  as  Ok  decreases  through 
Smax-  Consequently,  further  local  evolution  resembles  that  of  a  forward  parabolic  equation, 
serving  to  smooth  out  gradients  to  the  background,  sn  ~  scrit  on  the  fast  0(1)  time  scale. 
Moreover,  when  a  jump  07.  collapses,  then  that  grid  point  n, t  becomes  part  of  the  background, 
and  the  system  (7.20,  7.21)  is  reduced  to  a  (K  —  l)-dimensional  system  for  the  remaining 
jumps.  This  model  of  the  piecewise-in-time  evolution  of  system  (2.25)  is  supported  by  Figures 
10,  11  and  17a,  which  show  piecewise  smooth  dynamics  punctuated  by  the  collapses  of  jumps 
at  finite  times.  It  is  also  appropriate  to  note  that  the  reduced  model  (7.20)  gives  numerical 
results  that  are  indistinguishable  from  simulations  of  the  full  discrete  model  (7.13),  apart 
from  short  0(1)  transients  associated  with  jump  collapse  in  the  regime  scrjt  <  cr*  <  smax. 

The  finite-time  collapse  of  unstable  jumps  typically  occurs  one  jump  at  a  time,  (see  Fig¬ 
ure  10b).  In  a  further  simplification  of  the  reduced  model,  we  show  how  to  isolate  the 
evolution  of  a  single  jump,  and  show  that  the  simplification  leads  to  only  small  inaccuracies 
in  the  simulation,  and  the  possibility  of  increased  understanding  of  the  collapse  mechanism. 
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Fig.  19.  (left)  Comparison  of  the  piecewise  linear  flux  F(s),  (7.18),  predicted  by  (7.22)  after  the 
collapse  of  jump  a 6  with  the  actual  solution,  (right)  Comparison  of  the  decay  of  jump  predicted 
by  (7.22)  with  the  actual  evolution. 

Let  us  assume  that  in  the  coarsening  process,  there  is  a  separation  in  the  timescales  of  the 
successive  collapses  of  the  jumps  {a^}-  This  assumption  is  valid  if  there  is  a  separation  in  the 
values  of  52F{ok)  for  the  jumps  ak,  i.e.  if  —  S2F(oi)  ~  — <52-F(<72)  ~  —S2F(a3)  •  •  •  <C  —$2F(oj) 
then  Oj :  will  be  the  next  jump  to  collapse;  this  is  the  case  for  <r6  in  Figure  18b.  Then  during 
the  finite  time  that  the  jump  Uj  collapses  to  <jj  — >-  smax,  the  remaining  jumps  will  have 
changed  only  slightly  (see  Figure  19a).  Consequently,  if  we  neglect  the  slow  evolution  of  the 
other  jumps  sk,  then  system  (7.20)  reduces  to  a  single  first-order  ODE  for  <jj  while  all  of  the 
other  Ok  are  held  constant, 

d(Jj  =  Aj  ~  BjF{<jj) 

dt  Ax2 

where  the  constants  Aj,  Bj  are  given  by 

A  _  F(aj+ 1)  ,  F(aj- i)  B  _  nj+i  —  rij-i 

J  nj+l~n3  n3~n3~  1  ^  (nJ+l  ~nj)(n3  ~n3-lY 

Equation  (7.22)  may  be  integrated  starting  from  the  initial  size  of  slope  Sj  at  time  TK,  when 
the  solution  has  K  jumps,  to  determine  time  when  collapse  has  occurred,  Oj  =  smax, 

leaving  a  solution  with  K  —  1  jump  discontinuities.  If  further  we  assume  that  all  of  the  jumps 
are  equally  spaced,  with  (n,j+ 1  —  rij)  Ax  —  1/K  and  the  jumps  neighboring  Oj  are  very  near 
equilibrium,  F(<jj- 1)  =  F{oj+\)  —  T(scrit),  then  (7.22)  becomes 


d(J3 

dt 


(sait)  -  F(oj))  ~ 


2NK  sin  a  sin  <fi 


(7.24) 


where  the  second  approximation  results  from  the  asymptotics  of  F(s)  for  large  s,  (2.15). 
This  equation  has  the  approximate  solution 


dj(t )  ~  \JaNK  sin  a  sin  <fi  (TK_ i  —  t ),  TK  <  t  <  TK_ i,  (7-25) 

where  Tk_i  is  the  collapse  time.  From  Figure  19b,  we  see  that  (7.25)  qualitatively  captures 
the  nature  of  the  finite  time  collapse  of  the  jump.  Figure  19  illustrates  the  use  of  the  piecewise 
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linear  flux  approximation  (7.18)  and  equation  (7.22)  to  calculate  the  collapse  of  the  <r6  jump 
starting  from  the  initial  conditions  given  in  Figure  18. 

This  simplified  model  (7.22)  can  be  expected  to  approximately  describe  the  dynamics  in  some 
intermediate  regime,  1  <C  K  <C  N.  As  long  as  K  N,  the  piecewise  linear  approximation 
(7.18)  will  describe  the  flux,  but  as  K  -»  1,  (7.22)  can  not  hold,  because  there  will  be  strong 
coupling  between  jumps  due  to  the  boundary  condition  constraint  (7.15).  In  practice,  we 
have  observed  that  the  long-term  dynamics  are  very  sensitive  to  the  spatial  coupling  of  the 
jumps.  A  scaling  law  is  observed  for  the  number  of  jumps  K  as  a  function  of  time  for  systems 
with  large  N .  Equation  (7.22)  does  not  capture  this,  but  the  jump  diffusion  model  (7.20, 
7.21)  does  reproduce  this  behavior  of  the  full  system  (2.25). 

If  we  use  the  results  of  the  reduced  jump-diffusion  model  (7.22,  7.25)  with  the  initial  values 
for  the  jump  sizes  obtained  from  the  numerical  simulations  from  Section  6,  <7j(TK_ i)  = 
0(N/K 3/2).  Then  we  obtain  that  the  transition  time  from  K  to  K  —  1  jumps  is 

ATk  =  0(N/K 4).  (7.26) 

Consequently,  the  cumulative  time  until  only  K  jumps  remain  is  given  by  the  summation 

Tk  =  E  ATk  ~  o(N )  £;  r4  =  o(n/k 3),  (7.27) 

k=N  k=N 


where  the  second  summation  can  be  expressed  exactly  in  terms  of  the  polygamma  function 
[1]  as  J2  k~3  =  —  ip^(K  +  l))/6  ~  K~3  / 3  +  0(Ar_3,  AT-4).  This  result  agrees  with 

the  estimate  K  —  0((N /f)1/3)  from  (6.3).  It  is  not  clear  how  to  derive  the  scaling  result  for 
Gj(TK- 1)  from  (7.22)  or  (7.20). 

In  conclusion,  we  have  shown  that  in  the  continuum  limit,  the  slow  timescale  for  evolution 
in  (2.25)  diverges  as  the  microscopic  discretization  lengthscale  vanishes,  Ax  —  N~l  — >  0. 
This  singular  behavior  is  a  consequence  of  the  asymptotic  form  of  the  non-monotone  flux 
function,  F(s)  for  s  — >-  oo.  Further  work  focusing  on  the  influence  of  different  forms  of  the 
flux  function  F(s)  is  being  pursued. 
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